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ABSTRACT 


In  the  standard  treatment  of  stability  in  chemical 
reaction  systems  by  normal  mode  analysis,  the  stability  of 
the  non-equilibrium  steady  states  is  usually  determined 
numerically.  The  recent  development  of  a  general  analytic 
method  for  linear  stability  analysis  has  allowed  the 
asymptotic  boundary  of  the  unstable  steady  states  in  para¬ 
meter  space  to  be  determined,  without  assignment  of  the 
rate  constants.  A  diagrammatic  counterpart  of  this  method 
is  used  to  pinpoint  the  interactions  which  may  cause 
instability,  and  also  retains  information  on  matrix  quali¬ 
tative  stability. 

The  subject  of  the  stability  analysis  in  this  thesis 
is  the  Oregonator  model  which  was  proposed  to  explain  the 
significant  mechanistic  features  of  the  Belousov-Zhabotinski 
reaction.  The  model  is  analyzed  with  the  reverse  reactions 
included.  A  comparison  is  made  with  the  irreversible 
model  to  determine  the  effect  of  reversibility  on  the 
region  of  instability  in  the  network.  A  variable  stoichio¬ 
metric  coefficient  is  also  contained  in  the  model,  and  the 
effect  of  this  parameter  on  the  boundary  of  the  unstable 
region  is  discussed  in  detail. 
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CHAPTER  ONE 


Background 

1.1.  Introduction 

In  the  past  decade,  the  behavior  of  chemical  reactions 

beyond  the  equilibrium  steady  state  has  been  investigated 

with  increasing  enthusiasm. ^  This  interest  arose  primarily 

due  to  the  discovery  of  biological  clocks  and  oscillating 

biochemical  systems.  Such  systems  include  oscillatory 

metabolic  reactions  related  to  oxidative  phosphorylation 

and  glycolysis.  Several  review  articles  on  the  subject 

of  oscillating  reactions  are  available. ^  ^ 

Other  areas  of  interest  involving  oscillating 

phenonema  are  hydrodynamics , ^ ^  economics, ^  and  ecology. ^ 

Ecologists,  who  studied  the  cyclic  fluctuations  in 

animal  populations,  promoted  the  use  of  qualitative 

9  10 

criteria  for  stability  of  matrices.  '  Such  a  criter- 

9 

ion  was  provided  by  Quirk  and  Rupert  as  a  result  of 
studies  done  on  economic  systems.  The  matrix,  which 
determines  the  stability  of  an  ecological  system,  describes 
the  predator-prey  interactions  between  populations  and  is 
called  the  community  or  interaction  matrix.  The  inter¬ 
action  matrix  used  by  ecologists  is  analogous  to  the 
matrix  formed  in  the  linear  stability  of  chemical  systems. 

In  the  matrices  associated  with  chemical  systems,  every 
off-diagonal  element  of  the  matrix  of  the  linearized 
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stability  problem  contains  a  parameter  which  also  appears 
in  a  diagonal  element,  as  a  result  of  the  reaction 
stoichiometry .  "*"~*"  This  restricts  the 

possible  matrix  element  variation  and  means  that  matrices 
occurring  in  chemical  networks  can  violate  the  criterion 
for  matrix  qualitative  stability,  and  yet  have  no  regions 
of.  instability.  An  extension  of  the  work  of  Quirk  and 
Rupert  to  the  chemical  stability  problem  has  been  done 
by  Clarke,"*"*"  utilizing  the  quantitative  method  of 
analysis  used  in  this  thesis.  Such  qualitative  results"*""*" 
from  studies  of  chemical  systems  should  apply  to  inter¬ 
action  matrices  in  ecology  and  economics  also,  since 
the  parameters  of  these  systems  have  similar  interrela¬ 
tionships  . 

Other  qualitative  criteria  for  'stability  in 

chemical  systems  exist  for  the  special  case  in  which 

the  reactions  of  the  network  obey  mass  action  kinetics. 

12 

One  such  criterion,  due  to  Shear,  '  guarantees  the  unique¬ 
ness  and  global  stability  of  equilibrium  steady  states 
in  closed  systems,  as  long  as  mass  action  kinetics  hold. 
Since  equilibrium  and  the  states  close  to  equilibrium 
cannot  oscillate,  they  are  not  the  ones  of  interest  in 

this  analysis.  Another  criterion  is  provided  by  Horn's 

13 

zero  deficiency  theorem,  which  gives  a  sufficient 
condition  for  qualitative  stability  of  any  reaction 
system  obeying  mass  action  kinetics.  Since  there  is  no 


• 
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restriction  to  the  states  close  to  thermodynamic  equi¬ 
librium,  this  criterion  is  more  powerful  than  Shear's 
result.  In  particular,  Horn's  theorem  greatly  restricts 
the  number  of  mass  action  systems  that  can  possibly  be 
unstable.  However,  it  cannot  guarantee  instability 
will  result  in  any  particular  system  which  violates  the 
conditions  of  the  theorem,  in  the  model  to  be  analyzed  here 
some  steps  -  like  many  experimentally  determined  rate 
expressions  -  do  not  obey  mass-action  kinetics,  and 
therefore  Horn's  theorem  does  not  apply.  The  necessity 
for  considering  non-mass  action  systems  in  the  non¬ 
equilibrium  region,  has  motivated  other  approaches  to 
the  stability  of  oscillating  chemical  systems,  which 
shall  now  be  mentioned  briefly. 

In  general,  the  method  of  stability  analysis  used 
to  study  oscillating  reactions  can  be  classified  as 
either  thermodynamic  or  kinetic  in  approach. ^  The 
stability  analysis  of  chemical  systems  has  been  approach¬ 
ed  by  Prigogine  et  al .  from  a  thermodynamic  viewpoint 

by  a  local  formulation  which  centers  around  the  second 

14  15 

law  of  thermodynamics.  '  In  nonlinear  processes,  however, 
a  thermodynamic  potential  function  is  available  close  to 
equilibrium  only,  where  the  system  must  be  stable. 

Far  from  equilibrium, where  oscillations  may  occur,  a 
sufficient  condition  for  stability  can  be  expressed  in 
terms  of  the  second  variation  of  entropy  production. 


' 


. 
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even  though  no  potential  exists.  This  so  called 
"Glandsdorf f-Prigogine  criterion"  is  not  used  in  practice 
to  test  stability  because  chemical  systems  are  defined 
kinetically,  while  the  Glandsdorf f-Prigogine  criterion 
can  only  be  used  if  the  independent  variables  are  the 
fluxes  and  affinities  of  irreversible  thermodynamics. 

The  introduction  of  thermodynamics  into  the  stability 
problem  increases  the  complexity  of  the  problem  and  for 
this  reason  chemical  reactions  have  not  been  readily 
analyzed  by  these  methods. 

Another  thermodynamic  attempt  to  formulate  chemical 

reaction  dynamics,  called  network  thermodynamics,  has 

1 6 

been  based  on  electrical  network  theory.  In  this 
approach  chemical  reactions  are  represented  by  electric 
circuits  in  which  the  forces  and  flows  obey  Kirchhoff's 

i 

Laws.  The  aim  was  to  allow  classical  network  theory 
to  be  applied  to  chemical  systems  through  a  linear  graph 
theory.  The  important  underlying  idea  is  that  it  is  the 
way  in  which  the  reactions  are  linked  (i.e. ,  the  way  the 
circuit  is  hooked  up)  that  determines  the  stability  of 
the  network.  Stability  analysis  within  this  framework 
leads  to  the  Glandsdorf f-Prigogine  criterion.  As  before, 
the  main  drawback  is  that  the  independent  variables 
are  the  chemical  potentials  instead  of  the  dynamic 
concentrations.  For  this  reason,  network  thermodynamics 
has  the  same  drawbacks  as  the  older  thermodynamic  approach 


■ 
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of  Prigogine. 

In  the  kinetic  approach  to  chemical  stability 

problems,  the  dynamic  equations  of  motion  are  written 

in  terms  of  the  chemical  concentrations  and  the  reaction 

velocities  explicitly.  The  usual  procedure  in  a  kinetic 

investigation  is  to  apply  a  so  called  linear  stability 

3 

or  normal  mode  analysis.  Thrs  is  the  only  method  which 

allows  instability  conditions  to  be  determined  directly  in 

terms  of  the  system  parameters.  However,  the  standard 

approach  has  difficulties  in  dealing  with  large  chemical 

systems.  As  a  result,  a  new  approach  to  the  linear 

stability  analysis  of  chemical  systems  was  developed 
17 

by  Clarke,  to  enable  the  study  of  realistic  chemical 
systems.  Such  improvements  in  the  theory  were  necessary 
before  quantitative  analyses,  such  as'  the  one  presented 
in  this  thesis, were  possible.  The  essential  background 
for  this  theory  is  given  in  chapter  2,  while  the  standard 
treatment  of  stability  by  normal  mode  analysis  is  briefly 
introduced  in  section  1.2. 

The  basic  idea  used  in  the  approach  of  this  analysis 
is  that  chemical  networks  have  certain  topological 
properties  related  to  the  dynamics  of  the  system.  The 
topology  of  the  network  is  important  because  it  plays 
an  essential  role  in  the  algebraic  relationships  between 
the  rate  constants  and  concentrations,  which  determine 
the  boundary  of  unstable  states  in  parameter  space. 


6. 


This  boundary  is  useful  for  testing  experimental  reaction 

mechanisms  and  rate  constants.  Also,  the  validity  of 

models  which  represent  complicated  chemical  reactions 

can  be  assessed  by  examination  of  the  network  topology. 

An  important  conjecture  in  this  regard  is  that  networks 

with  similar  topology,  stoichiometry  and  kinetics  have 

1 8 

similar  stability  properties.  This  has  been  substan¬ 
tiated  in  several  cases  and  should  simplify  the  analysis 
of  complex  networks. 

In  this  paper  a  five  step  model  of  the  Belousov- 
Zhabotinski  reaction,  which  exhibits  oscillations,  is 
analyzed  with  the  reverse  reactions  included.  The 
boundary  of  the  unstable  steady  states  is  found  analytic¬ 
ally.  The  analytic  method  used  here  has  a  diagrammatic 
19 

counterpart  which  can  be  used  to  visually  recognize 
the  unstable  terms  in  the  interaction  matrix.  The 
interaction  matrix  is  the  matrix  of  coefficients  of  the 
linear  system  of  equations  consisting  of  the  perturbed 
kinetic  equations  linearized  about  steady  state.  The 
principal  aim  of  the  analysis  is  to  establish  relation¬ 
ships  between  the  unstable  non-equilibrium  steady  states 
and  the  responsible  destabilizing  terms  in  the  interaction 
matrix.  A  variable  stoichiometric  coefficient  is 
retained  in  the  matrix,  and  the  change  in  the  boundary 
of  the  unstable  states  with  this  parameter  is  also  found. 
The  effect  of  the  reverse  reactions  on  the  boundary  of 
the  unstable  region  is  similarly  of  special  interest. 


. 
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1.2.  Kinetic  Stability 

The  dynamics  of  a  chemically  reacting  system  is 
governed  by  the  kinetic  equations.  These  equations  are 
in  general  a  set  of  coupled  nonlinear  partial  differen¬ 
tial  equations,  which  can  be  written  as 


+  V2ci,  i 


(1.1) 


1 


.  .  .  ,  n 


where  c(r,t)  is  the  concentration  vector  whose  components 
c^  are  the  time  and  space  dependent  concentrations  of  each 

species.  (c)  represents  the  reaction  kinetics  for  each 

reaction  creating  or  destroying  species  c^.  The  rate  of 

change  of  the  concentration  of  each  species  depends  on  the 

concentrations  of  the  other  species  in  the  system  and  the 

i 

rate  constants  Kj .  The  diffusion  of  species  c^  is  given  by 
the  second  term  in  equation  (1.1),  where  is  the  diffusion 
coefficient  given  by  Fick's  law.  Temperature  variations 
and  external  forces  of  the  type  encountered  in  hydrodynamics 
are  excluded.  The  system  of  equations  (1.1)  is  autonomous, 
and  as  such  the  rate  constants  are  not  considered  functions 
of  time. 

The  stability  of  the  solutions  of  the  system  (1.1) 

is  determined  from  the  stability  theory  of  nonlinear 

differential  equations.  A  particular  solution  of  the 

equations  (1.1)  can  be  written  as  the  set  of  functions 

u^(t),  i  =  1,  ...,  n.  This  solution  is  asymptotically 

20  . 

stable  in  the  sense  of  Liapounov,  if  all  other  solutions. 


-v; 


v^ (t) ,  coming  near  it  approach  it  asymptotically.  In 

more  rigorous  terms,  the  solution  u^  (t)  is  stable  if, 

given  e  >  0  and  tQ ,  there  is  n  =  ri(e,  t  )  such  that  any 

solution  v.  (t)  for  which  |u.  (t  )  -  v.  (t  )  I  <  n  satisfies 

| ui  (t)  -  vi (t)  |  <  e  for  t  >  t  .  If  u^(t)  is  stable, 

and  in  addition  |u.  (t)  -  (t)  |  ->  0  as  t  ->  °°,  then  u^  (t) 

20 

is  asymptotically  stable. 

Other  definitions  of  stability  which  apply  to  auto¬ 
nomous  systems  concern  the  parametric  trajectories  or 
orbits  of  (1.1),  which  are  the  curves  represented  by  the 
solutions  of  (1.1)  in  concentration  space.  Closed  orbits 
are  the  curves  represented  by  periodic  solutions  of 

(1.1) ,  and  correspond  to  physical  oscillations.  If 

u^ (t)  is  a  nonconstant  periodic  solution  of  (1.1)  then 
u^  (t  +  to)  =  u^(t)  for  all  t,  where  to  is  some  minimum 
positive  constant,  the  period.  Let  C  be  an  orbit  of 

(1.1) .  Then  C  is  orbitally  stable  if,  given  e  >  0, 
there  is  r)  >  0  such  that,  if  R  is  a  representative  point 
of  another  trajectory  within  a  distance  e  of  C  at  time 
x,  then  R  remains  within  a  distance  n  of  C  for  t  >  x. 

If  no  such  n  exists,  C  is  orbitally  unstable.  If  C  is 

orbitally  stable  and  in  addition,  the  distance  between 

R  and  C  tends  to  0  as  t  then  C  is  asymptotically 

20 

orbitally  stable. 

An  orbit  which  is  asymptotically  orbitally  stable 
is  a  stable  limit  cycle.  A  limit  cycle  is  an  isolated 


.  i 
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9. 

closed  trajectory,  and  is  necessary  for  sustained  oscilla¬ 
tions  which  are  independent  of  initial  conditions.  A 
limit  cycle  is  a  periodic  motion  in  a  nonlinear,  noncon¬ 
servative  system  which  attracts  the  trajectories  nearby, 
in  the  same  way  that  a  singular  point  such  as  a  stationary 
state  or  equilibrium  state  attracts  its  nearby  trajec¬ 
tories  if  stable.  Another  important  property  of  some  limit 
cycle  type  oscillations  is  that  they  can  be  self-exciting 
or  self-starting.  If  a  system,  which  is  initially  at 
rest,  can  evolve  spontaneously  to  a  stable  limit  cycle  due  to 
the  infinitesimal  fluctuations  in  the  system  alone,  then 
the  system  is  said  to  exhibit  "soft  excitation",  while 
the  other  possibility,  a  "hard  excitation",  requires  an 
initial  impulse,  greater  than  some  threshold  value.  In 
two  dimensions,  soft-excitation  occurs  when  an  unstable 
stationary  state  (focus)  is  enclosed  by  a  stable  limit 
cycle.  Unfortunately,  in  higher  than  two  dimensions, 
the  relationship  between  limit  cycles  and  the  nature  of 
stationary  states  has  not  been  found. 

Since  the  solutions  of  equations  (1.1)  cannot  be 
easily  found  in  general,  a  starting  point  for  the 
investigation  of  limit  cycle  type  oscillation  in  chemical 
reactions  is  the  stability  of  steady  (stationary)  states. 

The  steady  states  of  equations  (1.1)  are  defined  by 

dc . 

tst  =0/  i  =  1 ,  .  .  .  ,  n  (1.2) 


, 

• 

' 

10. 

and  the  spacially  homogeneous  steady  state  solutions  are 

the  set  of  concentrations,  represented  by  the  components 

of  the  vector  c  =  (c,  ,  c0  , ,  , . . ,  c  ),  which  satisfy 

^o  lo  2o  no  1 

the  algebraic  equations 

F. (c)  =  0  ,  i  =  1,  .  .  .  ,  n „  (1.3) 

l  ~ 

From  the  preceding  definitions,  the  steady  state  cq 
will  be  asymptotically  stable  if  all  trajectories  start¬ 
ing  sufficiently  close  to  it  tend  to  it  asymptotically 
as  t  ->  «. 

If  the  system  (1.1)  is  closed  to  the  flow  of  matter, 

then  the  concentrations  of  the  species  evolve  until  the 

rate  of  change  of  reagents  into  products  and  vice  versa 

is  equal.  The  final  steady  state  reached  represents  a 

21 

state  of  chemical  equilibrium.  Thermodynamic  and 
12 

kinetic  arguments  both  show  that  these  steady  states 
are  always  stable.  This  means  that  an  arbitrary  fluc¬ 
tuation  in  a  closed  system  cannot  cause  the  system  to 
evolve  away  from  the  equilibrium  steady  state  or  the 
steady  states  that  are  close  to  thermodynamic  equilibrium. 

In  an  open  system,  a  flow  of  matter  through  the 
system  can  occur.  This  means  that  one  can  distinguish 
two  different  types  of  species:  the  external  species 
A ^  which  flow  into  and  out  of  the  system,  and  the  internal 
species  which  are  contained  in  the  system  and  are 

9 

often  called  the  chemical  intermediates.  The  external 


■ 


■ 
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species  can  be  kept  at  a  constant  concentration  by 
adding  reagents  and  removing  products  continuously, 
whereas  the  internal  species  have  unrestricted  concentra 
tions,  which  evolve  according  to  the  kinetic  equations 
(1.1).  Since  the  external  species  can  be  fixed  at  some 
constant  concentration,  they  may  be  incorporated  into 
the  rate  constants,  so  that  they  do  not  appear  explicit¬ 
ly  in  (1.1).  The  new  set  of  apparent  rate  constants  k 
for  the  system  consists  of  the  true  rate  constants  K 
multiplied  by  the  constant  concentration  of  the  appro¬ 
priate  external  species.  The  steady  state  solutions 
of  (1.1)  are  now  found  by  solving  (1.3)  for  the  steady 
state  concentrations  of  the  internal  species,  c^q,  for 
some  particular  set  of  apparent  rate  constants.  For 

i 

every  choice  of  the  apparent  rate  constants  a  different 
solution  can  be  found.  At  one  particular  ratio  of 
the  apparent  rate  constants,  the  equilibrium 
steady  state  occurs,  and  corresponds  to  no  flow  through 
the  system.  All  other  solutions  will  correspond  to 
non-equilibrium  steady  states  where  a  net  flow  through 
the  system  is  necessary  for  maintaining  steady  state. 

In  this  way  the  thermodynamic  equilibrium  steady 
state  solution  is  extended  into  the  non-equilibrium 
region  of  parameter  space.  When  the  system  is  constrain 
ed  far  enough  from  equilibrium,  unstable  steady  states 
may  occur.  This  instability  arises  because  the  kinetic 


. 


12. 


equations  are  not  linear  and  occurs  at  some  finite 

distance  from  equilibrium,  beyond  the  so-called  thermo- 
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dynamic  branch  of  solutions.  These  unstable  steady 
states  are  of  interest  because  most  models  of  oscillat¬ 
ing  chemical  reactions  are  associated  with  unstable 
steady  states  in  open  systems,  although  an  unstable 
steady  state  is  not  mathematically  necessary  for  oscil¬ 
latory  phenonema. 

The  stability  of  the  non-equilibrium  steady  states 
can  be  studied  by  perturbing  the  system  (1.1)  about  the 
steady  state.  Perturbed  solutions  of  the  form 

c  (r , t)  =  c  +  C(r,t)  (1.4) 

represent  the  deviation  of  the  system  away  from  the 
spacially  homogeneous  steady  state  concentration  vector 
cq.  The  initial  perturbation  Z,  (r,o)  is  assumed  to  be 
small  in  comparison  with  cq.  If  the  system  is  stable, 
it  is  to  be  expected  that  the  perturbed  solutions  c^(r,t) 
remain  near  the  steady  state  solutions  c^  as  time  evolves. 
This  is  tantamount  to  saying  that  the  fluctuations 
C(r,t)  die  out,  (i.e.,  C(r,t)  0  as  t  -*  °°)  . 

r\J  n-J 

The  perturbed  solutions  (1.4)  are  substituted  into 
equation  (1.1)  and  a  Taylor  expansion  of  the  resulting 
equation  is  made,  which  retains  the  terms  to  first  order 
in  the  perturbation.  This  yields  the  variational  equations 


. 
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I 

j=l 


o 


+  D. 
1 


i  =  1 ,  .  .  .  ,  n 


(1.5) 


which  are  valid  in  a  neighborhood  of  the  steady  state  if 
the  eigenvalues  of  the  resulting  matrix  of  the  coeffi¬ 
cients  do  not  have  zero  real  parts. 

In  the  case  that  an  eigenvalue  with  a  zero  real 
part  occurs,  the  nonlinear  terms  of  (1.1)  must  be  includ¬ 
ed  before  stability  can  be  determined.  These  solutions 
correspond  to  particular  values  of  the  constrained 
parameters  {K^,  A^}  in  a  transformed  parameter  space 
{h,  j}.  The  new  parameters  {h,  j}  are  defined  later  in 
section  2.2.  In  any  case,  the  nonlinear  terms  decide  only 
whether  the  boundary  is  included  in  the  unstable  region 
or  not,  and  the  linearized  equations  are  always  valid 
up  to  the  boundary.  It  is  this  boundary  which  will  be 
determined  in  this  analysis,  for  the  model  which  is 
introduced  in  the  next  section. 


' 


. 
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1.3.  The  Belousov-Zhabotinski  Reaction 

In  chemistry,  one  of  the  historically  important 

oscillating  reactions  is  the  Belousov-Zhabotinski  reaction, 

the  cerium  catalyzed  oxidation  of  malonic  acid  by  bromate 

in  aqueous  sulphuric  acid.  This  reaction  has  steady 

state  solutions  to  the  kinetic  equations,  beyond  an 

instability  of  the  thermodynamic  branch  of  solutions. 

These  steady  states  if  unstable,  can  give  rise  to  the 

oscillation  of  the  concentration  of  certain  reagents  in 

time.  Indeed,  this  is  the  case  as  observed  experiment- 
23 

ally.  In  addition,  this  reaction  exhibits  spacial 
structures  in  the  form  of  travelling  concentration  bands, 
when  the  temporal  oscillations  are  coupled  with 
diffusion. ^ ^ 

i 

The  great  interest  in  the  Zhabotinski  reaction 

24 

caused  the  formation  of  a  simple  model,  the  Oregonator, 
to  explain  the  important  features  of  the  mechanism. 

This  model  consists  of  the  following  five  steps: 

(1)  A  +  Y  +  X 

(2)  X  +  Y  2P 

(3)  B  +  X  ->  2X  +  2Z  (1*6) 

(4)  2X  +  Q 

(5)  2Z  +  fY 

where  A  =  B  =  [PrC>3“] ,  P  =  [HOBr] ,  X  =  [HBr02] ,  Y  =  [Br“], 

Z  =  [Ce (IV) ]  and  Q  =  B  +  P. 


' 


. 


$ 
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In  the  oriqinal  paper  Z  =  [2Ce(IV)]/ 

and  the  coefficients  of  2  in  front  of  Z  in  steps  (3) 

and  (5)  are  left  out.  Several  chemical  reactions  sum 

to  form  step  (5) ,  and  the  stoichiometric  coefficient  f 

represents  experimental  factors  which  are  not  yet 

determined.  The  actual  value  is  thought  to  be  close  to 

25 

1,  but  Jwo  and  Noyes  found  experimentally  that  f  varied 

continuously  in  the  range  0-2. 

The  Oregonator  model  (1.6)  has  been  subjected  to 

analysis  by  numerical  integration  techniques . ^ ^ 

However,  the  analysis  was  complicated  by  the  stiffness 

of  the  Oregonator  kinetic  equations.  A  good  discussion 

of  the  difficulties  encountered  in  numerical  integration 

27 

of  stiff  differential  equations  is  given  by  Gelmas. 

The  model  was  later  extended  to  include  the  reverse 
2  8 

reactions.  Due  to  the  oscillatory  nature  of  the 

system,  existence  of  a  limit  cycle  was  proposed  and 

29 

later  proved  analytically. 

The  reversible  model  of  (1.6)  possesses  a  unique 
homogeneous  stationary  state  which  becomes  unstable 
beyond  a  bifurcation  point.  The  appearance  of  a  second 
(periodic)  solution  to  the  kinetic  equations  then  gives 
rise  to  stable  homogeneous  oscillations.  These  oscilla¬ 
tions  exhibit  limit  cycle  behavior,  and  excitable  dynamic 
30 

Turner  has  demonstrated  the  existence  of  a  hysteresis 
cycle  between  the  stable  stationary  and  oscillatory 


states.  An  analysis  of  a  more  complicated  model  of  the 
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Zhabotinski  reaction  has  been  done  by  Clarke  '  with 
the  reverse  reactions  omitted. 

In  the  present  analysis,  the  Oregonator  model  (1.6) 
with  the  reverse  reactions  is  analyzed,  with  the  variable 
stoichiometric  coefficient  f  included.  Even  with  such  a 
simple  model,  the  boundary  between  the  stable  and  unstable 
steady  states  in  the  network  is  greatly  complicated  by 
the  presence  of  the  reverse  reactions.  Although  these 
results  are  interesting  in  themselves,  generalizations 
which  increase  the  understanding  of  chemical  network 
stability  are  also  sought  in  this  analysis. 


. 
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CHAPTER  TWO 


Stability  Theory  for  Chemical  Networks 


2.1.  Introduction 

In  this  chapter,  an  analytic  method  developed  by 
Clarke  for  determining  the  boundary  of  the  unstable 
steady  states  in  parameter  space  is  explained.  This 
theory  is  then  applied  to  the  model  (1.6)  and  the  results 
are  given  in  Chapter  3.  The  effect  of  diffusion  is 
neglected  and  only  the  homogeneous  perturbation  from 
steady  state  is  considered.  In  this  case  the  variational 
equations  (1.5)  are  given  in  matrix  notation  as 


(2.1) 


where 


and  £  =  c  -  cQ  represents  a  small  fluctuation  away  from 
the  steady  state  solution  c  .  M  is  called  the  inter- 

u  ~T>  = 


action  matrix,  and  characterizes  the  linearized  stability 
problem  about  steady  state.  The  solutions  of  the  matrix 
equation  (2.1)  in  the  case  where  |  is  a  constant  n  x  n 


matrix  are  found  from  the  matrix 
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(2.2) 


17 


. 


' 
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which  is  a  principal  matrix  of  solutions  of  (2.1).  To 
understand  the  evolution  of  any  initial  5,  it  is  suffi- 

r>J 

cient  to  calculate  the  matrix  (2.2). 

If  g  is  a  constant  nonsingular  matrix,  the  change 
of  coordinates 

C  =  P  x  (2.3) 

^  ^ 

yields  the  equivalent  system  of  differential  equations 

x  =  P-1  M  P  x  (2.4) 

^ 

which  can  be  written  as 

x  =  A  x  (2.5) 

where  A  is  similar  to  It  follows  that  the  matrices 

A  and  M  have  the  same  characteristic  polynomial 

det(Iw  -  M)  =  det(Io)  -  A)  (2.6) 

H  r*>j  — 

and  so  the  eigenvalues  are  also  the  same.  It  is  now 
sufficient  to  discuss  the  solutions  of  the  system  (2.5), 
whose  fundamental  matrix  is 


where  one  can  consider  that  g  is  such  that  ^  is  in 
Jordan  canonical  form. 

VJhen  in  Jordan  canonical  form,  ^  is  a  block  diagonal 
matrix  where  each  block  corresponds  to  one  eigenvalue  w . , 


. 


■ 
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and  each  block  has  the  structure 
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(2.7) 


where  B..,  B.„,  .  ..,  B.  are  basic  Jordan  blocks  belong- 
—  —  D  ^  — 

ing  to  the  eigenvalue  ok  .  A  basic  Jordan  block  correspond¬ 
ing  to  is  a  matrix  with  the  structure 


(2.8) 


Defined  in  this  way,  the  dimension  of  the  matrix  Aj 
corresponds  to  the  multiplicity ,  y,  of  the  eigenvalue 
a) j  .  The  number  of  basic  Jordan  blocks  r^  associated 
with  an  eigenvalue  ok  is  equal  to  the  nullity  of  the 
eigenvalue  w ^ .  The  nullity,  v,  is  defined  as  the  number 
of  linearly  independent  eigenvectors  which  correspond 


to  an  eigenvalue  of  multiplicity  y. 

A  ’  t 

The  matrix  e=3  can  be  visualized  by  considering 

B  ’  #  t 

the  structure  of  each  block  e^O1  which  can  be  written 


as 


. 
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The  solutions  of  (2.5)  are  linear  combinations  of  the 
columns  of  e=^  and  the  solutions  of  (2.1)  are  obtained 
from  (2.3) . 


If  the  eigenvalues  are  distinct  then  A  reduces  to 


a  diagonal  matrix  (each  has  one  element,  w ^ )  and 
At 

e=  has  the  form 


(2.10) 


The  solutions  of  (2.1) 

n 


C  (t) 

/V 


=  I 

k=l 


can  now  be  written  as 


(2.11) 


<frv  =  Pv  e 


“kt 


where 
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is  an  arbitrary  constant  and  is  a  column  vector  of 

P.  p,  is  the  eigenvector  corresponding  to  the  eigen- 

value  aj,  ,  and  satisfies 
k 


M  Pk  -  w  p 

—  ^  IS. 


(2.12) 


If  the  eigenvalues  are  not  distinct,  then  each  of 
the  matrices  ^ ^  which  are  of  dimension  s,  say,  contri¬ 


bute  s  linearly  independent  solutions  of  the  form 

k-1 
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■  03  -i  t 

d  .  .  ,  =  e  J 

Z j+m-1 


m 

z  - — 

k=l  (k-1) ! 


Pj+m-k'  m  =  1'  3  (2-13) 


to  equation  (2.11),  which  no  longer  has  n  linearly 

independent  eigenvectors.  Each  of  the  solutions  in 

(2.13)  corresponds  to  a  linearly  independent  solution 

of  the  degenerate  eigenvalue  clk  ,  given  in  terms  of  the 

3  6 

generalized  eigenvectors ,  Pj+m-k'  sYstem*  Each 

generalized  eigenvector  satisfies  the  relation 

-  “j>m  Pj+m-k  -  °-  (2-14> 

In  particular,  if  m  =  1,  then  equation  (2.12)  is  obtained, 

and  p .  is  an  eigenvector  corresponding  to  o).  . 

~  D  3 

When  the  stability  of  these  solutions  is  examined, 
it  turns  out  that  only  the  sign  of  the  eigenvalue  is  of 
importance,  since  the  factor  ewj^  dominates  the  series 
in  t.  However,  if  is  not  distinct  and  has  a  zero 


' 


. 
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real  part,  then  the  series  in  t  in  equation  (2.13)  is 
important.  The  problem  now  reduces  to  an  examination 
of  the  eigenvalues  of  £4,  determined  by  the  equation 

det(Io)  -  M)  =0.  (2.15) 

This  determinant  has  the  associated  characteristic 
polynomial , 

n  _ . 

p(w)  =  det  (Irn—  ^)  =  \  a .  con  1,  (2.16) 

i=0  1 

and  consequently  the  search  for  the  eigenvalues  of  M 
is  equivalent  to  finding  the  roots  of  the  real  polynomial 
(2.16) . 

In  conclusion,  the  eigenvalues  m  determine  the 
stability  of  the  general  steady  state  cq.  If  all  the 

i 

eigenvalues  have  negative  real  parts,  the  solution  is 
asymptotically  stable.  If  one  of  the  eigenvalues  has  a 
positive  real  part,  the  solution  is  unstable.  If  no 
eigenvalue  has  a  positive  real  part,  but  some  eigen¬ 
value  has  a  zero  real  part,  nothing  can  be  said  since 

20 

the  nonlinear  terms  in  (1.1)  affect  the  stability. 

The  eigenvalues  oo  are  a  function  of  the  rate 
constants  and  the  concentrations  of  the  reactants  and 
possibly  the  products.  Solution  of  the  stability 
problem  with  the  rate  constants  as  the  parameters  is 
sometimes  possible,  but  in  general  the  algebra  is 
considerably  simplified  if  a  new  set  of  independent 


■ 
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parameters  is  chosen.  These  parameters  are  designated 
h  and  j  and  are  introduced  in  section  2.2.  The  h,  j 
parameters  were  specifically  chosen  so  that  the  inter¬ 
action  matrix  M  could  be  written  as  a  bilinear  form  in 

h  and  j . 

~  ~ 

The  next  step  in  the  analysis  is  to  find  a  criterion 

which  can  guarantee  the  sign  of  the  real  parts  of  the 

eigenvalues  of  M.  Such  a  criterion  is  provided  by  the 

37 

Routh-Hurwitz  conditions,  which  are  the  subject  of 

section  2.3.  The  Routh-Hurwitz  criterion  is  based  upon 

the  signs  of  the  Hurwitz  determinants  A^,  i  =  1,  . . . ,  m, 

whose  elements  are  the  coefficients  a.  in  the  character- 

i 

istic  polynomial  of  ^  given  by  equation  (2.16);  the 

t  h 

i,  j  element  in  the  i  determinant  is 

—  a2  j  -i '  ^  ^ r  -j  <  &  7  if  2  j  -i  <  0  then 

“2j-i  =  °'  <2-17) 

and  the  criterion  states  that  the  real  parts  of  the  roots 
of  p(oo)  are  negative  provided  the  Hurwitz  determinants  are 
all  positive  definite.  Since  the  eigenvalues  w  change 
continuously  with  the  parameters,  the  boundary  between 
the  stable  and  unstable  states  occurs  when  the  largest 
real  part  of  the  eigenvalues  of  ^  is  zero.  Therefore, 
the  boundary  of  the  unstable  states  occurs  when  one  of 
the  Hurwitz  inequalities 

>  0  £  =  1,  ...,  m  (2.18) 

is  first  violated,  (i.e.,  A^  =  0  for  some  Z . ) 


. 

' 


At  this  stage,  the  stability  analysis  has  been 

reduced  to  the  problem  of  finding  the  zeroes  of  the 

Hurwitz  determinants,  which  are  homogeneous  polynomials 

in  the  elements  of  M.  The  approximate  zeroes  of  the 

Hurwitz  determinants  are  found  by  applying  a  method  for 

finding  the  asymptotes  of  real  power  polynomial  surfaces 

3  8 

which  has  been  recently  developed  by  Clarke.  The 

method  is  a  generalization  of  Newton's  Polygon  which  was 

first  used  in  1676  to  find  the  roots  of  two  variable 
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polynomials  by  geometric  considerations.  The  applica¬ 
tion  of  this  theory  to  the  chemical  stability  problem 
yields  an  equation  for  the  asymptotic  boundary  of  the 
unstable  states  in  parameter  space  and  is  discussed  in 
section  2.4.  In  addition  to  the  analytic  theory  of 
section  2.4,  a  diagrammatic  approach  is  used  to  summariz 
the  connection  between  the  terms  in  ^  and  the  boundary 
of  instability.  These  aspects  of  the  theory  are  intro¬ 
duced  in  section  2.5. 

The  method  used  here  can  be  applied  to  interaction 
matrices  outside  the  field  of  chemical  kinetics.  In 
particular,  an  application  was  made  to  the  Rayleigh- 
Benard  instability  of  a  binary  mixture.  The  complica¬ 
tions  introduced  in  the  more  general  case  will  not  be 
considered  in  this  thesis,  however. 


' 
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2.2.  Chemical  Network  Parameters 


The  chemical  stability  problem  is  usually  set  up 
with  the  rate  constants  as  the  independent  parameters. 
However,  in  this  case  the  interaction  matrix  usually  has 
elements  which  are  irrational  polynomials  in  the  rate 
parameters.  This  complication  can  be  avoided  by  a 
judicious  choice  of  a  new  set  of  independent  parameters. 

Given  the  rate  constant  k  for  each  reaction,  the 
steady  state  concentrations  and  the  rate  of  each  reac¬ 
tion  are  determined,  and  can  be  used  as  an  alternate 
set  of  independent  parameters  (c  ,  v°)  which  corresponds 
to  (k,  e)  where  e  labels  the  possibly  multiple  steady 
states.  Any  steady  state  velocity  vector  v°  may  be 
represented  as  a  sum  of  contributions  from  various  indepen- 
dent  overall  reaction  processes  conserving  the  given  steady 
state.  These  processes  form  a  reference  set  of  velocity 
vectors  called  the  extreme  currents .  The  coefficients 
in  this  sum  are  called  the  currents ,  j^,  and  form  a  new 
set  of  independent  parameters  for  the  system,  which  must 
be  positive  to  represent  a  physical  reaction  process. 

The  new  set  of  parameters  can  also  be  denoted  by 
(h,  j)  where  h  has  components  which  are  the  reciprocal 
of  the  steady  state  concentrations  and  j  has  components 


26. 


which  are  the  currents  of  the  system.  As  will  be  shown 
below,  this  choice  of  independent  parameters  reduces  the 
interaction  matrix  to  a  bilinear  form  in  (h,  j)  and  then  the 
expansion  of  the  characteristic  polynomial  of  M  is  homo¬ 
geneous  in  the  parameters  of  the  system: 


n 


M  =  l  jk  M 
k=l  K  " 


(k) 


where 


m 


(k)  = 

'ij 


s.<k)  h,  , 
ID  D 


(2.19) 


and  the  summation  is  over  the  n  currents,  j^,  of  the 
system. 

(k) 

The  matrix  2  has  elements  which  are  some  pure 
numbers.  The  important  point  to  note  is  that  the 
stability  of  any  state  (k,  e)  is  identical  to  the  stab¬ 
ility  of  any  of  the  associated  states  (h,  j).  Therefore, 
finding  the  region  of  instability  in  (h,  j)  space  will 
determine  the  region  of  k-space  which  is  unstable.  The 
rest  of  this  section  is  devoted  to  deriving  the  current 
parameters  specifically  and  illustrating  the  theory 
with  relevant  examples. 

The  equations  of  motion  (1.1)  for  any  chemical 
system  can  be  written  in  matrix  notation  as 

dc 


(2.20) 


' 

< 
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where  v  is  a  constant  matrix  reflecting  stoichiometric 
effects  in  the  kinetic  equations,  and  v  =  (v-^, 
v  )  where  is  the  velocity  of  the  k^  "  reaction. 

If  there  are  n  independent  internal  species,  (i.e.,  the 
dependent  species  have  been  eliminated)  then  ^  is  a 
n  x  r  constant  matrix.  The  steady  state  condition 
becomes 


V  •  v  =  0  ,  (2.21) 

and  the  velocity  vectors  which  satisfy  this  equation 
will  be  denoted  by  v° .  Any  general  steady  state  of  the 
network  is  represented  by  a  vector  . 

Since  the  n  species  are  independent,  the  rows  of 
y  span  a  n-dimensional  subspace  S  which  is  orthogonal 
to  the  r-n  dimensional  complementary  subspace  Sj_,  in 
which  v°  lies.  All  this  means,  is  that  the  solutions 

/■v 

v°  must  be  orthogonal  to  the  subspace  in  which  v  lies. 

As  a  result,  the  steady  states  v°  can  be  expressed  in 
terms  of  any  basis  set  which  spans  Sj_: 


l 

k=l 


(2.22) 


where  the  are  some  arbitrary  constants.  The  minimal 
basis  set  for  Sj_  will  contain  r-n  basis  vectors. 

However,  there  is  a  directional  requirement  on  the 
velocity  of  a  reaction. 


9 


i  =  1  / 


•  •  •  9 


r 


(2.23) 


' 


-* 

which  is  necessary  if  the  reaction  is  to  occur  in  the 
direction  as  written.  A  component  of  v  which  is  negative 
has  no  physical  meaning:  it  is 

not  equivalent  to  the  reverse  of  a  reaction.  Take  for 
example  the  reaction 


A  +  X 


Y  +  Z. 


The  velocity  of  the  forward  reaction  is  K^AX  while  the 
velocity  of  the  reverse  reaction  is  KrYZ .  These  rates 
are  only  equal  at  equilibrium.  At  a  non-equilibrium 
steady  state  (open  system)  the  rates  are  not  equal,  and 
it  is  this  lack  of  symmetry  in  the  flow  of  matter  that 
complicates  the  dynamics  of  chemical  systems.  In  a 
reversible  network,  each  reaction  and  its  reverse  must 
be  treated  as  two  separate  reactions,  each  with  its  own 
non-negative  velocity  v^. 

The  conditions  (2.23)  on  the  velocity  components 
place  a  restriction  on  the  choice  of  the  basis  set  of 
vectors  i,  because  v°  must  now  be  expressed  as  a  posi- 
tive  linear  combination  of  i: 


(2. 


Each  of  the  conditions  (2.23)  determines  a  half space 
which  passes  through  the  origin.  The  solutions  of  (2.21) 
will  have  to  lie  in  a  subspace  of  the  convex  polyhedral 


' 


. 

■ 
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cone  defined  by  the  intersection  of  the  halfspaces 
(2.23).  With  this  restriction,  the  correct  set  of  basi 
vectors  (2.24)  necessary  to  span  the  solution  space  of 


v  are  found,  although  this  basis  is  not  necessarily 
the  minimal  one. 

When  the  steady  state  velocity  vector  is  expressed 
as  in  (2.24),  the  are  a  positive,  unrestricted  set 
of  parameters  which  satisfy  the  steady  state  condition 
(2.21)  subject  to  (2.23).  The  j,  in  (2.24)  are  there¬ 
fore  the  required  current  parameters.  The  associated 

set  of  basis  vectors  {i,  }  of  Si  are  called  the  extreme 

~k 

currents  and  are  a  reference  set  of  currents  for  the 
29 

network. 

The  complete  set  of  extreme  currents  defines  a  set 
of  halflines  which  extend  from  the  origin  in  the  direc¬ 
tion  given  by  each  of  the  vectors  i^.  This  set  of  half 
lines  determines  a  set  of  halfspaces  whose  intersection 
is  a  convex  cone,  denoted  by  C  .  The  set  of  half lines 
which  form  the  edges  of  this  cone  are  said  to  span  the 


cone  and  together  they  define  the  cone's  frame. 


40 


The 


points  which  lie  on  the  surface  or  interior  of  C  are 

c 

the  physically  acceptable  solutions  of  (2.21),  v°  e  Cc, 

and  satisfy  (2.24).  C  is  called  the  system's  steady 

c 

29  . 

state  current  cone,  and  lies  m  the  r-n  dimensional 


subspace  Sj_. 


.. 
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The  intersection  of  the  current  cone 
unit  hyperplane, 


C  with  the 
c 


v  •  (1 ,  1 ,  . . . ,  1)  1  ,  (2.25) 

yields  a  convex  polytope,  called  the  current  poly  tope 

II  .  All  steady  states  on  a  line  through  the  origin  in 
c 

C  are  mapped  into  a  single  point  in  II  ,  and  the  extreme 
c  c 

points  of  the  current  polytope  are  the  extreme  currents. 

Any  general  steady  state  velocity  vector  v°  of  the  net- 

rsJ 

work  is  given  by 

^  =  X  \  jk  ik  '  (2'26) 

k 

where 


k 


1  , 


0  , 


and  A  is  a  scaling  factor. 

The  current  polytope  is  a  cross-section  of  the 
current  cone  and  lies  in  a  subspace  of  one  less  dimen¬ 
sion  (r-n-1) .  The  important  thing  to  note  is  that  all 

steady  states  in  C  on  a  line  through  the  origin  have 

c 

the  same  stability  by  virtue  of  (2.19),  and  so  need  not 
be  distinguished.  The  significance  of  the  current 
polytope  is  that  it  can  be  used  to  decompose  the  full 


■ 
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stability  problem  into  a  sum  of  smaller  ones,  by  a 
procedure  called  simplicial  decomposition . ^  This  method 
is  typically  used  for  reducing  the  dimension  of  the  stability 
problem  when  the  complete  set  of  extreme  currents  is  larger 
than  the  minimal  basis  set,  as  often  occurs  in  practice. 

For  example,  take  the  current  polytope  consisting  of  four 
points  in  the  plane  as  follows: 


j2 


where  each  of  the  points  j  ^  to  corresponds  to  an 
extreme  current  of  the  network.  Either  of  the  diagonal 
lines  joining  or  j2^4  can  ke  use<^  to  form  two 

simplicies  (triangles)  consisting  of  three  points.  The 
stability  problem  can  now  be  solved  for  any  point  in 


. 
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the  plane  by  using  the  current  parameters  determined 
by  the  simplex  in  which  it  lies.  This  reduces  the 
number  of  current  parameters  needed  in  the  interaction 
matrix  by  one. 

In  order  to  illustrate  the  concept  of  an  extreme 

current,  two  simple  examples  will  be  discussed  next. 

Then  an  extreme  current  of  the  Oregonator  model  (1.6) 

will  be  derived  as  an  example  of  those  which  are  found 

in  the  results  of  Chapter  3.  As  the  first  example, 

consider  a  reaction  network  consisting  of  a  stepwise 

transformation  of  matter: 

K1  K  K  K4 

A - - — ►X-j- - =— -*'X2 - - — *-X3 - - — ►  B .  (2.27) 

The  X^  are  the  internal  species  and  the  A,  B  are  external 
species.  Kj^,  ...,  K4  are  the  rate  constants  of  the 
reaction  steps  with  corresponding  velocity  v^,  ...,  v^ . 

It  can  be  seen  immediately  that  only  one  steady  state 
exists  in  this  system,  and  the  requirement  on  the  rates 
v^  is  that 

V1  =  v2  =  v3  =  v4  • 

The  matrix  of  the  net  stoichiometric  coefficients  for 
(2.27)  is 


' 
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and  the  equations  of  motion  are 


dx 


where  v  can  be  written  explicitly  in  terms  of  its 

components  as  =  K-^A,  v2  =  K2X1 ,  =  K^X2 ,  =  K^X^ 

th 

and  X^  is  the  concentration  of  the  i  species.  The 
steady  state  solutions  are  the  velocity  vectors  v°  which 
satisfy 

v  •  v°  ~  0  ,  v?  0  .  (2.28) 

=  ~  l 

The  rows  of  y  span  a  three  dimensional  subspace  of  the 
four  dimensions  of  v.  Therefore,  solutions  of  (2.28) 

'■'W 

are  the  vectors  v°  which  lie  in  the  one  dimensional 

subspace  orthogonal  to  y.  The  current  cone  C  for  this 

system  is  a  single  vector  which  extends  in  the  direction 

i^  =  (1,  1,  1,  1)  from  the  origin,  where  i^  is  the 

extreme  current  for  the  network  (2.27).  The  intersection 

of  Cc  with  the  unit  hyperplane  given  by  v°*  (1,  1,  1,  1)  =  1, 

yields  the  current  polytope  II  for  the  system,  which  in 

c 

this  case  is  a  single  point,  representing  the  extreme 
current  i..  This  is  the  simplest  convex  set  consisting 
of  one  extreme  point  only.  The  general  steady  state 
solution  is  a  constant  multiple  of  i.-^, 


.  . 
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where  is  an  arbitrary  constant. 

As  a  second  example,  consider  the  reaction  sequence 


Once  again  the  internal  species  are  denoted  by  X^,  while 
the  external  species  are  A,  B,  C.  In  this  case  two 
alternate  extreme  reaction  sequences  are  seen  which  can 
satisfy  the  steady  state  requirement: 

v,  =  v_  =  vc  =  =  v_ 

1  3  5  6  7 

v2  =  v4  =  v5  =  v6  =  v?  . 

Here  the  subscript  on  corresponds  to  the  marked 
in  (2.29).  These  sequences  correspond  to  two  extreme 
currents : 


11  =  a,  o,  i,  o,  i,  if  i) 

12  =  (0 ,  1,  0,  1,  1,  1,  1) 

which  are  obtained  by  solving  (2.28)  with  the  net 
stoichiometric  matrix 


v 


1  0 
/  0  1 

I  0  0 

y  o  o 

\  0  0 


-10  0 
0-10 
11-1 
0  0  1 

0  0 


0  1 


' 
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and  where  v  =  (v^  v  ,  v3  ,  v  ,  v5,  v  ,  v?). 

In  this  example,  the  null  space  of  (2.28)  is  a 

2-dimensional  subspace,  and  the  extreme  currents  i,  and 

i_2  determine  two  halflines  which  span  the  current  cone 

C  which  lies  in  a  plane: 
c 


(2.30) 


Once  again,  the  intersection  of  the  current  cone  with 

the  unit  hyperplane  yields  the  current  polytope  II  for 

c 

the  system,  which  in  this  case  is  a  line  segment: 


■ 


\ 


r 
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(2.31) 


whose  endpoints  are  the  extreme  currents  and  . 

Any  general  steady  state  solution  for  the  system 
(2.29)  lies  on  this  line  segment  according  to 

Z°  =  x  A{  +  i'2  ^  J  0  ’ 

where  X  >,  0,  is  some  scaling  factor  and  =  1. 

In  other  words  (2.25)  has  the  effect  of  projecting  the 

points  which  lie  in  the  plane  (2.30)  into  the  line  in 

(2.31).  This  is  valid  because  the  steady  states  lying 

on  a  line  through  the  origin  in  C  have  the  same  stability. 

c 

The  final  example  will  be  drawn  from  the  Oregonator 
model  (1.6).  The  extreme  current  derived  here  occurs 
in  both  the  reversible  and  irreversible  models,  and  is 


\ 


one  of  the  important  sources  of  instability  in  the 
analysis  presented  in  Chapter  3.  For  this  reason,  the 
treatment  will  be  detailed.  From  the  Oregonator  reac¬ 
tions  in  (1.6),  the  equations  of  motion  are 

X  =  k1Y  -  k2XY  +  k3X  -  2k4X2 

Y  =  -k1Y  -  k2XY  +  fk5Z  (2. 

Z  =  2k3X  -  2k5Z 


where  the  external  reactants  A  and  B  are  included  in  the 
rate  constants  k.,  so  that  the  network  can  be  considered 

l 

an  open  system  depending  on  the  three  dynamic  concentra¬ 
tions  X,  Y  and  Z.  The  k's  are  written  explicitly  as 
ki  =  K^A,  k2  =  I<2 ,  k3  =  K  B,  k^  =  I<4  and  k^  =  . 

In  the  usual  matrix  notation,  these  equations 
become 

dc 


where  c  =  (X,  Y,  Z)  is  the  concentration  vector,  v  is 
the  net  stoichiometric  matrix. 


v 


1  -2 

0  0 

2  0 


/ 


and  v  is  the  velocity  vector  whose  components  are 


' 


1 


V 


• 
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2 

v,  =  k.Y,  v~  =  k~XY,  v_.  =  k_X,  v„  =  k.X  ,  vc  = 

1  12  2  3  3  4  4  5 

Every  steady  state  satisfies 

v  •  v°  =  0  ,  v?  >.  0  . 

=  ~  l 

The  rows  of  y  span  a  three  dimensional  subspace  S  which 
is  orthogonal  to  the  complementary  2-dimensional  sub¬ 
space  Sj_  in  which  v°  lies.  It  follows  from  the  preced¬ 
ing  discussion  that  the  current  cone  C  is  2-dimensional 

c 

and  that  the  current  polytope  II  is  one  dimensional. 

c 

The  net  result  is  that  only  two  extreme  currents  can 
exist  in  this  case  and  they  are  the  extreme  points  of 
the  line  segment  IIc. 

Two  features  of  this  model  are  different  from  the 
usual  case.  First  the  stoichiometric  parameter  f  in  reaction 
(5)  is  included  in  v,  and  secondly  the  reactions  (1.6)  are 
not  linearly  independent.  In  this  example,  only  the  case  f>l 
is  considered,  and  further  details  are  left  until  Chapter 
3.  The  singular  nature  of  y  is  overcome  by  utilizing 
the  condition  for  normalization  of  the  velocity  in  an 
extreme  current, 


in  addition  to  (2.28).  The  condition  (2.34)  provides 
an  extra  equation 


(2.33) 


(2.28) 


39. 


(1,  1,  1,  1,  1)  •  v°  =  1 


and  can  be  combined  with  (2.28)  to  yield 


v '  v°  —  b 


(2.35) 


where 


1-1  1-2  0 


-1-1  0  0  f 


0  0  2  0  -2 


11111 


and  the  transpose  of  b'  is  the  vector  (0,  0,  0,  1). 

One  may  show  that  the  extreme  solutions  of  (2.35)  are 
found  by  deleting  all  possible  combinations  of  columns 

i 

(reactions)  of  y'  until  no  steady  state  remains.  When 
the  particular  columns  1,  2,  3,  and  5  are  chosen,  equa¬ 
tion  (2.35)  becomes 


-1  -1  0 


1 


0  0  2 


o 


0 


1  1 


Solving  these  relations  for  v°  yields  the  extreme  current 


o 


f-1  f+1 


1,0,1),  f  £  1 


2 


2 


r 


(2.36) 
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and  the  scaling  factor  l/f+2,  where  i'  =  l/f+2  i  from 
(2.26)  . 

The  extreme  current  represents  a  steady  state 
of  the  system  in  which 


vi  = 
o 

V2  = 
o 

V3  = 
o 

= 


f-1 


=  k,Y 
1  o 


mi  =  k,xoY 

2  2  o  o 


1  =  k3Xo 


1  =  kcZ 

5  o 


for  the  current  =  1  M/sec.  The  mapping  between  the 

rate  constants  and  the  steady  state  concentrations  X  , 

Y  ,  Z  for  i,  can  be  written  from  the  equations  above  as 
o  o  ~A  ^ 


 (f-1) 


 (f+1) 


k-,  =  — ,  k0  =  ,  k~  =  — ,  ,  kc  =  —  . 


2Y 


o 


2X  Y 
o  o 


1_ 

X 


1_ 

z 


Substituting  these  relations  into  equations  (2.32)  now 
yields  the  following  set  of  dynamic  equations: 


At  this  stage  it  can  be  seen  that  the  extreme  current 
represents  the  full  nonlinear  kinetic  equations  for 
the  reactions  (1) ,  (2) ,  (3)  and  (5)  of  the  Oregonator 


model  (1.6) . 


. 


' 
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The  interaction  matrix  M,  can  be  written  for  this 

=A 

extreme  current  as  follows.  Consider  an  arbitrary 

« 

perturbation  £  =  (AX,  AY,  AZ)  =  (X-XQ,  Y-Yq ,  Z-ZQ)  about 
the  steady  state  (Xq,  Yq ,  Zq)  and  neglect  the  nonlinear 
terms,  then  the  equations 

AX  =  AY/Yq  -  AY/Yq  -  AX/Xq  +  AX/XC 

AY  -  -  fir)  ay/yo  -  (r r )  AY/Yo  -  (rr)  ax/xo  +  f  AZ/Zo 

AZ  =  2  AX/X  -  2  AZ/Z 
o  o 


are  obtained.  If  new  variables  are  defined  by  h  =  1/X  . 

J  x  o' 

h  =  1/Y  ,  h  =  1/Z  ,  then  these  equations  are  equivalent 
V  o  z  o 

to  equation  (2.1)  where  is  given  by  (2.19)  for 
is  written  explicitly  as  follows: 


f  >  1  .  (2.37) 


Each  extreme  current  i^  of  the  system  has  an  associated 
matrix  which  appears  in  the  interaction  matrix  ^  as 
a  weighted  sum,  j^M^,  as  stated  in  equation  (2.19). 

For  the  irreversible  model  (1.6)  two  more  extreme 
currents  were  similarly  found: 


(0,  f,  1,  ,  1) 


1+f  ,  1) 
2 


f  4  1 


(2.38) 


(f,  0,  1, 


all  f. 


■ 
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In  each  region  of  f,  the  elements  of  the  interaction 
matrix  M  are  a  sum  of  the  elements  occurring  in  each 
matrix  corresponding  to  an  extreme  current  ,  as 
follows : 


jA 

SA 

+ 

jC 

Sc  ' 

f 

1 

^  B 

=B 

+ 

=c  ' 

f 

< 

1. 

Note  that  i.  =  .i  at  f  =  1  as  required  by  continuity  of 
the  solution  of  f  =  1. 

The  above  procedure  for  obtaining  the  interaction 

matrix  from  the  net  stoichiometric  matrix  is  completely 

analytic.  In  the  case  of  complex  reaction  networks,  a 

computer  program  exists  which  calculates  all  the  extreme 
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currents,  by  the  procedure  given  above.  The  representa¬ 
tion  of  the  general  steady  state  of  the  network  as  a 
positive  linear  combination  of  the  extreme  currents, 
guarantees  that  the  correct  overall  stoichiometry  is 
obtained  for  the  net  chemical  change  of  the  external 
species.  To  illustrate  this  important  result,  the 
irreversible  Oregonator  model  will  be  used  as  an  example 
of  the  difficulties  which  may  arise  in  finding  the 
correct  overall  reaction  in  a  network. 

The  original  claim  that  the  Oregonator  generates 

the  unique  stoichiometry  fA  +  2B  2P  +  Q  for  the  net 

43 

reaction,  has  been  corrected  by  Noyes.  In  his 


. 


discussion,  the  correct  overall  reaction  is  given  as  a 
positive  linear  combination  of  any  two  of  the  following 
steps,  except  0c  or  0^ 

0  (f-l)A  +  2B  +  (l+f)P  ,  f  >  1 

a 

0b  2B  +  2fP  +  (l-f)Q  ,  f  <  1 

0  2f  A  +  2B  ■>  (1+f )  Q  ,  all  f 

0,  B  ->  P 

d 

C>  fA  +  2B  +  fP  +  Q 

e 

0f  2 (f-1) A  +  2B  +  2P  +  (f-l)Q,  f  >  1 

in  the  system's  oscillatory  region. 

The  steps  0^  to  0^  represent  some  of  the  possible 
overall  reactions  which  can  occur  when  the  network  (1.6) 
is  at  steady  state.  It  can  be  shown,  using  the  extreme 
currents  for  this  network  that  only  the  steps  0  to  O 

ci  C 

are  necessary,  and  any  general  steady  state  of  the  net¬ 
work  is  a  positive  linear  combination  of  either  0^  and 
0  or  0,  and  0  .  Since  any  general  steady  state  in  the 
network  is  a  positive  linear  combination  of  the  extreme 
currents,  the  general  overall  reaction  can  be  written 
as  a  positive  linear  combination  of  the  extreme  overall 
reactions.  The  extreme  overall  reactions  are  obtained 
directly  from  the  extreme  currents  which  represent  a 
particular  linear  combination  of  the  reactions  in  the 
network.  For  example,  2i.^  =  (f-1,  f+1,  2,  0,  2)  from 
equation  (2.36)  represents  a  possible  steady  state  of 


' 

' 


' 
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the  network  consisting  of  a  sum  of  reactions  (1),  (2), 

(3)  and  (5)  as  follows: 


(f-1)  (Rl)  +  (f+1)  (R2)  +  2  (R3 )  +  2  (R5) 

where  i.  is  multiplied  by  two  to  remove  fractions.  When 
the  reactions  in  (1.6)  are  substituted  into  this  sum, 
the  extreme  overall  reaction  given  in  step  C>a  above  is 
obtained.  Similarly,  when  the  extreme  currents  i^  and 
i_  given  in  (2.38)  are  considered,  the  extreme  overall 
reactions  found,  can  be  identified  with  the  stoichio¬ 
metries  C>k  and  C>c  above,  respectively.  Since  any 
general  steady  state  of  the  network  is  a  positive  linear 

combination  of  i„  and  i_  for  f  >  1,  step  0,-  above  can 

~A  ~C  ^  f 

be  expressed  in  terms  of  and  d.  as  follows 


f  >  i. 


(2.39 


The  step  0^  which  is  valid  in  both  regions  of  f  represents 
the  following  linear  combination  of  the  extreme  currents: 


f  >  1 


4- 


t 


f 


< 


1 


+ 


+ 


/ 


f  =  1 


Finally,  the  step  0^  represents  the  extreme  overall 

reaction  for  the  extreme  currents  i.  or  in  at  f  =  1. 

~A  ~B 

(Recall  i  =  L  at  f  =  1.) 

~A  ~B 


. 

In  summary,  the  extreme  overall  reactions  are  0^, 
0^  and  C>c .  Any  general  steady  state  of  the  network 
(1.6)  has  an  overall  reaction  which  is  some  positive 
linear  combination  of  the  extreme  overall  reactions;  in 
particular  for  each  region  of  f : 


^A°A  +  ^c°c  '  3a  +  jc  =  1. 


3B°B  +  Jc°c  '  3B  +  3C  =  1, 


f  i-  1 
f  4  1 


where  is  the  current  corresponding  to  the  extreme 
current  i^.  Since  A  =  B  =  [BrO^  ] ,  B  can  be  eliminated 
from  the  overall  reactions  represented  by  these  relations 
as  follows: 


( ^A  +  2jc)A  jA  P  +  jcQ  '  f  >  1 
{ 2  j  B  +  2(l+f)  jc}A  ->•  2f  jB  P  +  {  (1-f )  jB-  +  (1+f)  jc}  Q,  f  «  1 

The  ratio  of  the  currents  for  the  network  (1.6)  can  be 
determined  from  these  equations  if  it  is  known  whether 
f  is  greater  or  less  than  one,  for  example: 

j  :j  is  in  the  ratio  P:Q  f  ^1. 

Any  particular  steady  state  of  the  network  is 
represented  by  a  certain  ratio  of  currents,  and  has  a 
corresponding  overall  reaction  which  is  a  sum  of  extreme 
overall  reactions.  The  extreme  currents  have  determined 


the  overall  reaction  stoichiometry  without  any  assumptions 


. 


. 


- 

' 


being  necessary. 
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A  O 

Using  the  relative  rates  of  the  reactions,  Noyes^ 
concluded  0^  and  0^  were  not  permissible  overall  reac¬ 
tions  in  the  systems  oscillatory  region.  From  the 
analysis  of  the  Oregonator  model  (1.6)  in  Chapter  3,  it 

will  be  shown  that  i  is  a  stable  extreme  current  and  so 

~c 

0  would  represent  the  extreme  overall  reaction  of  a 
stable  steady  state  for  any  value  of  the  rate  constants. 

Step  O^,  on  the  other  hand,  contains  the  extreme  current 
i-A  which  is  the  source  of  instability  for  f  >.  1 .  The 
stability  of  the  steady  states  represented  by  the  extreme 
overall  reaction  0^  will  therefore  depend  on  the  instability 
conditions  defining  the  boundary  of  the  unstable  region, 
which  will  give  the  necessary  condition  on  the  ratio  j  / j 

I  V 

for  which  the  linear  combination  j_i,  +  j  i  is  an  unstable 

JA~A  J c^c 

steady  state  for  any  particular  set  of  rate  constants. 


2.3.  Application  of  the  Stability  Criterion 
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Once  the  interaction  matrix  M  has  been  derived, 
the  characteristic  polynomial  (2.16)  is  known.  The 
network  is  stable  if  all  perturbations  from  the  steady 
state  decay,  or  alternately,  the  eigenvalues  of  M  have 
negative  real  parts.  This  can  only  occur  if  all  the 
Routh-Hurwitz  conditions  are  satisfied.  The  theorem  of 
Routh-Hurwitz  will  now  be  given,  as  well  as  the  explicit 
form  of  the  Hurwitz  determinants,  A^,  defined  in  equa¬ 
tion  (2.17) . 


37 

The  criterion  of  Routh-Hurwitz  states  that  all 
the  roots  of  the  real  polynomial  p  (to)  =  otocon  +  a-j  oon  ^  + 
...  +  an  (aQ  ^  0)  have  negative  real  parts  if  and  only 
if  the  inequalities 

{a  A  >  0 ,  n  odd 
o  n 

A  >  0,  n  even 
n 


all  hold.  The  A's  can  be  written  explicitly  if  a  >  0 

o 

as  follows: 


o 

A 

“l 

a3 

o 

A 

“l 

a3 

a5 

** 

• 

• 

• 

o 
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al 
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CXt-  . 
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In  the  case  aQ  >  0,  an  alternate  set  of  necessary 

and  sufficient  conditions,  the  Lienard-Chipart  stability 

37 

criterion,  is  available.  This  theorem  states  that 

for  all  the  roots  of  the  real  polynomial  p(u>)  =  aQ  l on  + 

n_ 

+  ...  +  an(a0  >0)  to  have  negative  real  parts, 
the  necessary  and  sufficient  conditions  can  take  any  of 
the  following  forms: 


i) 

a 

n 

> 

0, 

a  0 
n-2 

> 

0 

/  •  •  •  / 

A 

i — 1 

< 

0,  A3 

>  0 ,  ... 

ii) 

an 

> 

0, 

a  0 

n-2 

> 

0 

/  •  •  •  7 

a2  > 

0,  A4 

>  0 ,  ... 

iii) 

an 

> 

0, 

a  . 
n-1 

>0 

/ 

a  o>0 

n-3 

/  •  •  ®  / 

A1>0, 

A3>0 ,  .  .  . 

iv) 

a 

n 

> 

0, 

a  , 
n-1 

>0 

r 

a  _>0 

n-3 

/  •••  7 

a2>o. 

• 

• 

• 

o 

A 

<3 

The  advantage  in  using  one  of  these  sets  of  condi¬ 
tions  is  that  the  analysis  of  a  higher  order  Hurwitz 
determinant  is  avoided.  (In  this  analysis,  A^.). 

When  the  Lienard-Chipart  conditions  are  applied  to  the 

characteristic  polynomial  (2.16)  for  M,  a  set  of  stability 

inequalities,  f  (p)  >  0,  p  =  1,  ...,  k,  are  obtained 

from  (2.40)  which  are  polynomials  in  the  parameters  of  the 

system,  p  =  (h,j).  Network  instability  occurs  in  regions  of 

parameter  space  where  at  least  one  of  the  inequalities 

f  (p)  >  0  is  violated.  Therefore,  the  region  of  parameter 
P  ~ 

space  where  any  of  the  polynomials  f  (p)  is  negative  is 

p  ~ 

part  of  the  system's  unstable  region,  and  the  various 
portions  of  its  boundary  are  given  by  f  (p)  =  0  for  some 

p  ~ 

p.  A  method  for  finding  the  approximate  solutions  to 


the  equations  f  (p)  =  0  has  been  developed  by  Clarke 

p  ~ 

and  is  the  subject  of  the  next  section.  Since  the 
method  applies  to  every  stability  inequality,  the 
subscript  p  will  be  dropped  from  f  (p)  in  the  future 

P  r>*/ 
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2.4. 


Conditions  for  Instability 


The  unstable  region  of  the  network,  denoted  by 
f (p)  <  0  in  the  previous  section,  is  bounded  in  parameter 
space  by  the  surface  defined  by  f  (p)  =  0.  The  necessary 
approximation  to  this  surface  is  obtained  as  a  set  of 
asymptotes  in  the  theory  which  follows.  First,  however, 
an  example  from  the  Oregonator  model  (1.6)  will  be  used 
to  show  the  motivation  behind  the  development  of  such  a 
theory . 

We  shall  consider  the  polynomial  obtained  for 
from  the  interaction  matrix  M  given  in  (2.37).  For 
convenience,  consider  the  case  f  =  1.1,  then  the 
characteristic  equation  for  becomes 

(to  +  5.0  x  10  3x  y  0  \ 

1.05  x  to  +  l.ly  -l.lz  I  (2. 

-2x  0  to  +  2z  / 


where  new  variables  x=jh  ,  y  =  j  h  ,  z=jh  have 

A  X  y  Z 

been  defined.  Expanding  this  determinant  yields  the 
characteristic  polynomial 


P  (co) 


to3  +  to3  {5.0  x  10  3x  +  l.ly  + 
co{2.2yz  +  O.lxz  -  0.995xy}  + 
to  +  to  +  toa2  +  to  . 


2z  }  + 
to°  (0.21 


xyz) 


The  second  Hurwitz  determinant,  is  the  polynomial 


N 


formed  by  the  product  cuo^ 


since  aQ  =  1 ,  and  yields 


f-^(p)  =  ^2  (P)  =  4.4yz2  +  0.2xz2  +  2.42  y2z  -  1.81  xy2 

(2.42) 

-1.98xyz  +  5.0  x  10-2x2z  -  4.98  x  10  2x2y. 

In  order  to  illustrate  the  solutions  to  A?  (p)  =  0, 

^  'N/ 

the  projection  on  the  (log  x,  log  z)  plane  is  taken,  by 
setting  y  =  1  in  (2.42): 

f o  (p)  =  4.4z2  +  0.2xz2  +  2.42z  -  1.09x  -  1.98  xz 

Z  <N/ 


+  5  x  10  2x2z 


-  4.98  x  10  2x2 


(2.43) 


When  the  zeroes  of  this  polynomial  are  plotted  in  the  log^ 

parameter  space  tt  =  (log  x,  log  z),  the  curve  in  figure 

2.1  results.  The  asymptotes  to  this  curve  for  large  tt 

are  also  shown  in  figure  2.1. 

* 

If  a  vector  y  =  (y, ,  y0)  is  defined  corresponding 
to  each  term  of  f 2 (p) ,  where  y^  is  the  power  to  which  x 
is  raised  in  f 0  (p) ,  while  y0  is  the  power  to  which  z  is 
raised,  then  the  points  labelled  y,  to  yn  in  figure  2.1 
are  obtained.  For  example,  the  third  term  in  (2.43)  is 
2.42  x°z1  which  has  the  associated  vector  y^  =  (0,  1).  These 
seven  points  define  a  polygon  of  minimum  area,  which  is 
outlined  in  figure  2.1. 

The  interesting  feature  of  this  polygon  is  that  each 

asymptote  to  the  curve  in  figure  (2.1)  is  perpendicular 

*Do  not  confuse  the  vector  y  with  the  current  parameter 
y  =  1/Y  previously  set  equal  to  one. 
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Figure  2.1 

Plot  of  A0  (p)  for  the  extreme  current  i-,  at  f  =  1.1  of  the 

/  rv/ 

irreversible  model,  projected  on  the  tt  =  (log  x,  log  z) 
plane . 
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to  one  of  the  polygon's  edges.  The  second  point  to  note 
is  that  each  of  these  edges  contains  two  points  of  the 
polygon  which  correspond  to  terms  of  f „  (p)  with  coef- 
ficients  of  opposite  sign.  For  example,  from  figure  2.1 
the  points  associated  with  negative  terms  in  f 9  (p)  are 
yA,  Yc*  Ynr  an<^  the  asymptotes  are  perpendicular  to  the 
edges  containing  the  points  y9,  y.  and  y r,  y7. 

In  general,  these  observations  hold  for  any  real 
power  polynomial  and  are  the  basis  of  the  method  for 
finding  the  approximate  surface  of  f (p)  =  0  in  parameter 
space.  In  order  to  find  an  equation  which  will  yield 
the  asymptotes  in  figure  2.1,  the  ideas  of  the  preced¬ 
ing  example  will  now  be  put  on  a  firm  basis. 

Consider  any  stability  inequality  f (p)  >  0  which 
can  be  written  as  a  polynomial  in  the  parameter  set. 


I  0  <  p±  <  »,  i  =  1, 
as  follows 


f  (P) 

rv 


m 


I 

j=l 


ci  n  pkYkj 
3  k=l  K 


n}  , 


(2.44) 


where  the  polynomial  contains  m  terms  with  real  coefficients 

c-,  and  p  e  Rn  is  a  n-tuple  of  real  numbers.  The  itk 
D  ~ 

term  in  (2.44)  is 

c.p.^li  p„^2i  p_^3i,  ...,  p  ^ni 
1*1  ^2  ^3  ^n 


and  has  the  exponent  vector 


. 


•  * 


•  -  ^ 
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yni} 


(2.45) 


By  choosing  any  base  a,  a  e  R,  a  >  1,  for  the  logarithmic 
function,  the  term  p^yk  maY  be  written  as 


a 


lo9a  Pk>  . 


The  definitions 


Y  =  {y .  |  i  =  1 ,  . . . ,  n} 

X l  =  (l°9a  PX/  loga  p2,  log^  pn)  (2.46) 

t  h 

applied  to  the  l  term  of  (2.44)  now  gives 
it  •  y . 

c .  a~  , 

l 

where  it  and  y.  can  be  interpreted  as  points  in  the 

~  ~i 

Euclidean  space  En,  with  the  usual  definitions  holding 
for  the  inner  product  and  norm.  The  polynomial  (2.44) 
written  in  terms  of  the  transformed  parameters  tt  becomes 

g  (  tt  )  =  \  c(y)a~  y  (2.47) 

yeY  ~ 


where  c (y)  is  the  sum  of  the  coefficients  of  all  the 

/N/ 

terms  in  f (p)  with  the  exponent  vector  y: 


c  (y) 

rv 


I 


[yi  = 


y] 


(yi} 


(2.48) 


The  set  of  exponent  vectors  (y|c(y)  ^  0}  for  any 
g(u)  is  a  finite  point  set  Y  in  En .  This  set  of  points 
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determines  a  convex  polytope  called  the  exponent  poly  tope 
II,  which  is  defined  as  the  convex  hull  of  Y.  The  convex 
hull  of  a  finite  point  set  Y  is  the  smallest  convex  set 
containing  Y ,  defined  by  the  set  of  all  convex  linear 
combinations  of  y, ,  . . . ,  y  ,  which  is  the  set  of  line 
segments40 


m 


I  •  >  0 , 

i=l  1  1 


m 

1 

i=l 


X  . 

l 


=  1} 


(2.49) 


Returning  to  the  previous  example,  it  can  be  seen 
that  the  seven  points  in  the  polygon  of  figure  2.1 
correspond  to  the  exponent  vectors  of  f 0  (p)  given  in 
(2.43).  In  this  case 

Y2={(0,2),  (1,2),  (0,1),  (1,0),  (1,1),  (2,1),  (2,0)} 

I 

and  the  convex  hull  of  this  finite  point  set  is  the 
polygon  which  is  plotted  in  figure  2.1.  The  convex  hull 
of  Y2  consists  of  all  the  points  on  the  boundary  of  the 
polygon  plus  its  interior.  The  line  segment  joining  any 
two  points  of  Y2  can  be  written  as 

{ 7T 1 7T  =  Xy9  +  (1-X )  y,  ,  0  <  X  <  1}  (2.50) 

from  the  2-dimensional  case  of  (2.49).  The  line  segment 
(2.50)  is  visualized  below  on  the  line  AB  with  the  vector 
tt  restricted  to  the  line  segment  joining  y.  and  y9 

^  rs/  _L  'V/  Z 


44 


■ 
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When  X  is  unrestricted  (all  real  X) ,  the  set  of  points 

ir  in  (2.50)  is  the  line  passing  through  y.  and  y0. 

2 

In  E  ,  the  equation  for  a  line  can  be  written  as 

b  =  '  (2.51) 

where  b,  c-^  and  c2  are  given  constants.  Generalizing  to 
E  defines  a  plane;  b  =  +  c27r2  +  and  finally 

in  En  the  set  of  points  tt  satisfying 
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b  =  C,TT..  +  C~TT0  +  ...  +  C  TT  =  C  •  IT  (2.52) 

1122  n  n  ~  ~ 

are  said  to  define  a  hyperplane  for  given  values  of  b 
and  c.  Note  the  hyperplane  (2.52)  is  orthogonal  to  c. 

In  any  convex  polytope  II,  there  will  be  certain 
points  tt  which  do  not  lie  between  any  other  two  points 
of  the  set.  From  the  diagram  above  this  means  that  tt 
cannot  be  expressed  as  a  linear  combination  of  any  y. , 
y0  in  the  set  Y  with  0  <  X  <  1.  Or  in  other  words,  tt 
cannot  be  expressed  as  a  proper  convex  linear  combination 
of  any  other  two  points  in  the  set  IT.  In  this  case  tt 
is  called  an  extreme  point  of  the  set.  In  the  example 
shown  in  figure  2.1,  all  the  points  labelled  yn  to  y_, 
except  Yc,  are  extreme  points  of  the  exponent  polytope 
(hexagon) . 

From  the  discussion  above,  it  follows  that  an 
extreme  point  is  a  special  type  of  boundary  point  of  a 
convex  polytope.  In  particular,  the  extreme  points 
determine  the  surface  structure  of  II,  since  the  convex 
hull  of  the  set  of  extreme  points  in  Y  is  II.  In  general 
II  is  a  solid  embedded  in  a  n-dimensional  parameter  space, 
tt  =  (loga  p±,  loga  p2,  ...,  loga  pR)  .  The  surface 
structure  of  this  solid  can  be  described  as  sets  of  0, 

1,  2,  ...,  n-1  dimensional  faces,  which  correspond  to 

extreme  points,  lines  (edges),  planes  and  hyperplanes. 

The  significance  of  the  facial  structure  of  n  is 
that  it  can  be  used  to  form  a  series  of  approximations 


to  the  surf  ace  where  g  (tt)  =0.  In  other  words,  the 

polynomials  formed  from  the  terms  corresponding  to  the 

« 

points  in  Y  which  lie  on  the  faces  of  II,  denoted  by  g  ( tt  )  , 

r  ~ 

are  a  limiting  form  for  g  (n) .  To  show  this  result,  the 
generalization  of  (2.50)  to  the  case  of  a  half line  which 
starts  at  a  base  point  ttq  and  extends  to  infinity  in  the 
direction  tt  is  needed.  Such  a  halfline  is  called  a  ray, 
and  is  defined  as  the  set  of  points 

A  (tt  ,  TT)  =  {tt  (X  )  I  TT  (X  )  =  TT  +  X  TT ,  X  >  0  ,  TT(X)sEn}  (2.53 

1  7  /v  V  7  I  r>u  7  ' 

where  the  substitutions  yn  =  tt  ,  y~  -  yn  =  tt  have  been 

x-l  ~o  ~2  £  1 

made  in  (2.50).  Consider  any  point  tt(X)  on  a  ray  which 
is  perpendicular  to  the  hyperplane  passing  through  a 
face  F,  say,  of  II.  Assume  also  that  tt(X)  is  chosen  so 
that  gF(TTQ)  7^  0,  then  the  following  theorem  holds: 

3  8 

Theorem  1  (Clarke) 

lim  (3L (X)  )/g  (tt  (X)  )  =  1 

X->°° 

If  the  dimension  of  the  face  used  is  k,  then  g^  (tt )  is 
called  a  k-face  approximation.  In  the  chemical  stability 
problem,  the  approximation  usually  used  is  called  the 
edge  approximation,  because  the  points  y  on  the  1-dimen¬ 
sional  faces  of  II  are  used  to  define  gF  (tt)  . 

Theorem  1  allows  g  (tt)  to  be  replaced  by  the  poly¬ 
nomials  g1-,(TT)  in  the  limit  X+°°,  and  means  that  the  surface 
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g  (tt )  =  0  has  an  approximation  given  by  g^  (tt )  =0.  In 
the  case  of  an  edge  approximation,  where  no  interior 
points  occur  on  the  edge,  (2.47)  yields 

gF  (tt)  =  c±  a~*~l  +  c 2  a~’^2  (2.54) 

where  y. ,  y0  are  the  extreme  points  of  the  edge  and  cn , 
q-2  are  the  corresponding  coefficients.  If  the  zeros  of 
the  polynomial  g(ir)  =  0  are  to  be  approximated  by  gTI,(7r)  =  0, 
then  c^  and  c^  in  (2.54)  must  be  of  opposite  sign,  since 
the  exponents  are  always  >,  0.  Roughly,  what  one  is 
doing  is  finding  the  zeros  of  the  polynomial  g(Tr)  =  0 
by  examining  the  approximate  polynomials  gF  (tt)  =  0  formed 
from  the  points  of  Y  which  lie  on  an  edge  of  II,  where 
the  function  g(ir)  =  0  changes  sign. 

It  follows  that  using  higher  dimensional  faces  of 

II  for  gr,('n)  would  include  a  greater  number  of  terms  from 
r  *** 

g  ( tt )  in  the  approximation,  hence  improving  the  asymptotic 
surface  obtained  for  g(Tr)  =  0.  However,  the  edge  approx¬ 
imation  is  good  for  large  I  I  tt  I  I  ,  when  tt  is  not  almost 
orthogonal  to  a  higher  dimensional  face  of  II  (such  as  a 
plane) .  In  chemical  systems  where  the  concentrations 
of  intermediates  may  vary  over  large  orders  of  magnitude, 
the  condition  I  I  tt  I  I  large  is  usually  met  and  so  the  error 
incurred  for  small  I  |  tt  I  |  should  be  a  physically  insignif- 
leant  region  of  the  available  parameter  space. 


■ 
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The  equation  for  the  asymptotes  of  the  surface 
g(Tr)  =  0,  can  now  be  found  with  the  aid  of  two  more 
theorems  which  clarify  the  form  of  the  ray, 
which  is  necessary  if  AOr  ,tt)  is  to  be  an  asymptote  of 
g  (tt  )  =  o. 


3  8 

Theorem  2  (Clarke) 

Consider  the  direction  vector  tt,  and  the  highest  dimen¬ 
sional  face  F  of  II  lying  farthest  in  the  direction  tt. 
Choose  any  base  point  ^tto,  then  the  ray  A(jrQ,TT)  which 
extends  from  in  the  direction  tt  is  perpendicular  to 
F.  If  tt  is  such  that  g„  (tt  )  ^  0,  then  the  ray  A  (tt  ,tt  ) 
is  not  an  asymptote  of  the  surface  g  (tt )  =  0. 


This  theorem  is  actually  saying  that  asymptotic  rays  of 
the  surface  g  (tt  )  =  0  must  be  orthogonal  to  a  k-face  for 
k  >  0.  The  following  theorem  gives  a  sufficient  but  not 
necessary  condition  for  A(TTQ,n)  to  be  an  asymptote,  and 
can  be  used  to  find  a  base  point  ttq  for  an  asymptote  of 
g  (tt)  =  0. 


' 
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Theorem  3  (Clarke) ^ 

Consider  the  direction  vector  if,  and  the  highest  dimensional 
face  F  of  II  lying  farthest  in  the  direction  tt  .  Given 
the  ray  A(tto,tt)  which  is  perpendicular  to  F,  and  any  jto 
satisfying  gF(TrQ)  =  0  such  that  gF  (tt)  changes  sign  at  ^tto; 
then  the  ray  A(tto,tt)  is  an  asymptote  of  the  surface 
g  (tt)  =  0  as 


This  theorem  allows  an  asymptotic  ray  to  the  surface 
g(ir)  =  0  to  be  chosen  so  that  the  base  point  tt  lies  on 
the  edge  F ,  say,  and  satisfies  gF(T^Q)  =  0.  Hence  from 


(2.54) 

,  \  +  TT„*y,  .  -  tt  *y~  A 

g„  (tt  )  =  c  a^o  ^1  +  c  a^o  ^2  =  '0 

F  r->  O 


(2.55) 


where  y.  ,  y0  are  the  extreme  points  of  the  edge  F,  and 
c+ ,  c  are  the  positive  and  negative  coefficients, 
respectively.  Because  c+  >  0  and  c  <0,  gF(TTQ)  changes 
sign  at  ^tto  as  required  by  theorem  3  and  it  follows  from 

9f{£o}  =  0  that 

=  I  c+/c"  I 


or  taking  the  loga  °f  each  side 

=  loga|c+/c‘| 


(2.56) 
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The  component  of  perpendicular  to  the  edge  vector 
(y^y-^)  does  not  affect  this  equation.  If  any  ray 


A  (JL0  /  tt )  passing  through  is  chosen  so  that  tt  is  perpen¬ 


dicular  to  (y0-y,  )  also,  then  any  point  tt  =  tt  +  Att, 

A  >  0 ,  on  this  ray  will  also  satisfy  (2.56).  Therefore 
the  equation 


(2.57) 


holds . 


This  is  the  equation  for  a  hyperplane  (compare  with 
(2.52))  which  is  perpendicular  to  the  edge  of  II  joining 
the  extreme  points  y-,  and  y9.  The  hyperplane  (2.57) 
determines  an  asymptote  to  the  surface  in  II  where  g  (tt)  =  0, 
and  has  an  associated  halfspace  defined  by 


(2.58) 


which  gives  the  region  where  g  (tt)  is  negative.  Each  extreme 

point  of  II  is  called  a  positive  or  negative  vertex 

according  to  whether  the  corresponding  coefficient  is 

positive  or  negative.  Two  adjacent  extreme  points  are 

connected  by  an  edge  of  II.  The  intersection 

of  the  halfspaces  (2.58)  defined  by  each  edge 

of  a  negative  vertex  determines  a  polyhedral  cone  C  , 

which  gives  the  asymptotic  region  of  instability  due  to 

the  vertex.  In  any  exponent  polytope  II,  the  approximate 

k 

unstable  region  will  be  a  union  of  the  cones  C  defined 


Vv 
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by  each  negative  vertex.  The  total  approximate  region 

of  instability  in  a  chemical  network  will  be  a  union  of 

k 

each  set  of  cones  C  occurring  for  each  stability  inequality 
f(p)  >  0. 

Returning  to  the  example  of  figure  2.1,  the  equa¬ 
tions  for  the  asymptotes  can  be  determined  from  (2.57), 
for  the  two  edges  joining  y~,  y.  and  y a ,  y in  II.  This 

rJ  J  /v  4r  U  /V  / 

yields 


TT  * 

% 

-^3> 

=  TT  * 

/V 

df 

-1) 

=  log 

0.35 

TT  * 

/N/ 

~7 

=  TT  * 

(0, 

-1) 

=  log 

rH 

• 

O 

where  tt  =  (log  x,  log  z)  ,  and  the  conditions  for  instability 
(2.58)  become 

log  x  >  log  z  +  0.35 
log  z  <  1 

where  in  the  original  parameter  space  x  =  1/XQ,  z  =  1/Z0 
in  terms  of  the  steady  state  concentrations  XQ ,  Zq. 

It  should  be  noted  that  if  a  point  occurs  on  the 
edge  joining  adjacent  positive  and  negative  vertices, 
then  the  point  would  have  to  be  included  in  (2.55).  Let 
the  point  on  the  edge,  denoted  by  the  exponent  vector 
y^,  correspond  to  a  positive  coefficient  c->,  then 

(tt  )  =0  means  that 
^o 


c 


a 


+ 

c  a 


TT 

~o 


TT  *  V 

+  c^a  ~o  ~3 


from  (2.55).  If  y^  is  a  midpoint  of  the  ^2  ~  Xl  e<^9e  / 


■ 
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then 


£  *  (y o  -  y.i  )  _  + 


c  I  a  ~'o  vi2  XI '  -  c '  +  c^a  ~o  'X3 


£„•  (y.t  -  y  1 ) 


7T  •  (v  ““V  ) 

Let  a  Aj 3  XI'  =  £,  then  since  y0  -  y-,  =  2  (y_  -  yn  ) 

r+s  A  /~>S  -L  Jj 

it  follows  that: 


-  2  + 
c  r,  +c^?  +  c  =  0, 


and  the  equation  for  the  asymptotic  hyperplane  to  the 
edge  y0  -  y,  becomes,  in  parallel  to  (2.57), 

A.  r->  -L 

tt  *  (y0  -  y-. )  =  2  log  £ 

where  £  is  the  positive  root  of  the  quadratic  equation 
above : 

C  =  't”c3±  y/^3  +  4c_c+  }  l/2c_  . 


; 
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2.5.  Diagrammatic  Aids 

The  method  developed  in  the  last  section  allows  the 
boundary  of  the  unstable  states  to  be  determined  analytic¬ 
ally.  This  computation  can  be  further  simplified  if  the 
destabilizing  terms  can  be  recognized  by  inspection  of 
the  interaction  matrix  M.  In  order  to  accomplish  this, 

M  is  given  in  terms  of  a  collection  of  diagrams.  It  is 
the  purpose  of  this  section  to  show  that  these  diagrams 
are  equivalent  to  the  terms  of  M  and  that  certain  of  the 
cycle  diagrams  are  associated  with  instability. 

The  starting  point  for  the  diagrammatic  approach  is 

19 

the  reaction  diagram.'  The  reaction  diagram  contains 
each  chemical  reaction  in  the  network  together  with  the 
stoichiometry  and  kinetics,  and  corresponds  to  the  full 

i 

nonlinear  kinetic  equations.  For  example,  the  Oregonator 
model  (1.6)  has  the  following  reaction  diagram: 


' 
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Each  arrow  represents  a  reaction,  and  is  labelled  (1)  to 
(5)  in  a  parallel  fashion  to  the  reaction  steps  in  (1.6). 
Reactions  which  produce  or  consume  more  than  one  species 
are  denoted  by  a  split  arrow.  The  internal  species  are 
labelled  X,  Y  and  Z.  The  external  species  A,  B,  P  and 
Q  do  not  appear  in  the  reaction  diagram  as  they  can  be 
included  in  the  rate  constants.  An  arrow  which  points 
away  from  a  species  represents  a  reaction  which  consumes 
that  species,  while  incoming  arrows  represent  reactions 
which  produce  the  species.  The  number  of  feathers  and 
barbs  is  also  important.  A  standard  convention  for 
arrows  when  facing  from  the  tail  to  the  tip  is  that  the 
number  of  feathers  on  the  left  hand  side  of  the  shaft 
is  denoted  by  f  .  The  value  of  f  gives  the  kinetic 
order  of  the  reaction,  while  the  total'  number  of  feathers 
on  a  shaft  gives  the  stoichiometry.  The  number  of  barbs 
corresponds  to  the  stoichiometry  of  the  product,  and 
there  is  no  left-right  distinction.  As  an  example, 
consider  reaction  (3)  in  the  above  reaction  diagram, 
which  corresponds  to  the  reaction,  X  ->  2Z  +  2X. 

From  the  reaction  diagram  above  two  rough  cycles 
can  be  seen.  These  consist  of  R(l)  +  R(3)  +  R ( 5 )  and 
R(2)  +  R(3)  +  R(5).  Thus  one  may  suspect  that  two 
extreme  currents  exist  corresponding  to  these  combinations 
of  reactions.  In  the  actual  solution,  R(4)  will  also 
be  needed  to  conserve  the  balance  of  mass,  since  the 


' 

' 
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production  of  each  species  in  an  extreme  current  is  zero 
at  steady  state. 

As  a  specific  example,  take  the  extreme  current 
given  in  (2.36)  for  f  ^  1,  as 

,f-l  f+1  i  n  -i  \ 
i,A  2  '  2  '  1'  0'  1>  * 

This  vector  expresses  the  steady  state  solution  in  terms 
of  the  components  of  v°  which  satisfy 

(2/f-l)v°  =  (2/l+f)v°  =  v°  =  v°  .  (2.59) 

If  the  velocity  of  steady  state  in  reaction  (5)  is  set 
equal  to  j^,  v°  =  then  (2.59)  says  that  v^  =  (f-l)/2  j^. 

Therefore,  the  velocity  in  each  reaction  at  steady  state 
can  be  written  as  a  scaled  multiple  of  the  current  j  . 

This  can  be  written  diagrammatically  as  follows 


Z 


. 
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where  the  appropriate  reactions  from  the  reaction  diagram 
have  been  chosen.  This  diagram  is  called  an  extreme 
current  diagram .  It  is  interpreted  in  the  same  way  as 
the  reaction  diagram  given  previously.  When  a  label 
such  as  (f-l)/2  appears  beside  an  arrow,  the  entire 
reaction  is  multiplied  by  that  factor.  Hence  reaction 
(1)  is  (f-1  )/2  Y  ->  (f-l)/2  X. 

The  next  step  is  to  reduce  the  extreme  current 
diagram  to  the  interaction  diagram ,  which  corresponds  to 
the  linearized  perturbed  equations.  However,  the  inter¬ 
action  diagram  for  the  extreme  current  i  will  simply 
be  given  and  its  term  by  term  equivalence  to  the  inter¬ 
action  matrix  M  ,  given  in  (2.37),  shown.  The  interaction 
diagram  for  the  extreme  current  f  >  1 ,  given  above  is 


■ 

■ h  ■' 
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The  arrows  in  this  diagram  link  the  different  states 
of  a  reaction.  A  solid  arrow  links  a  reactant  to  a 
product,  and  corresponds  to  a  positive  matrix  element. 

A  dashed  arrow,  on  the  other  hand,  links  different  react¬ 
ants,  and  corresponds  to  a  negative  matrix  element. 

Solid  arrows  are  called  activators  and  dashed  arrows 
are  called  inhibitors . 

The  arrows  form  cycles  which  are  called  k-cycles, 
where  k  is  the  number  of  arrows  in  the  cycle.  For 
example,  the  interaction  diagram  above  contains  three 
cycles : 

1-cycle 


2-cycle 


3-cycle 


The  factors  f,  (f+l)/2  multiply  the  value  of  the  arrow. 
Each  reactant  has  an  associated  parameter  x  -  j.h  , 

A  X 

y  =  j  h  ,  or  z  =  j,h  ,  where  h  ,  h  ,  h,  are  the  reciprocal 
steady  state  concentrations  and  jA  is  the  current.  The 


' 
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value  of  an  arrow  is  given  by  the  reactant  at  the  arrow's 
tail,  multiplied  by  the  number  of  barbs  and  also  any 
factor  which  multiplies  the  arrow.  For  example,  the 
cycles  above  have  the  values  2x,  (f+l)/2  xy  and  -2fxyz, 

respectively.  The  sign  of  each  term  is  determined  by 
the  product  of  the  signs  of  each  arrow  in  the  cycle,  where 
dashed  arrows  are  negative  and  solid  arrows  are  positive. 

The  arrows  point  towards  or  away  from  a  second  type 
of  diagrammatic  element  called  a  stabilizer .  The  stabilizers 
are  denoted  in  the  interaction  diagram  above  as  (f+l)/2  x, 
fy  and  2z,  and  are  represented  diagrammatically  by  a  dot. 

The  stabilizers  contain  the  contribution  to  the  diagonal 
matrix  elements  coming  from  reactions  that  consume  a 
species,  and  are  negative  terms  in' Besides  the 
stabilizers,  the  only  other  diagrams  which  contribute  to 
iru^  are  the  1-cycles.  The  stabilizer  and  1-cycles  which 
occur  at  a  species  can  be  combined  to  form  a  modified 
coefficient,  while  the  diagram  which  survives  cancella¬ 
tion  is  kept  in  the  interaction  diagram  (i.e.,  a  stabil¬ 
izer  if  the  coefficient  is  positive,  or  a  1-cycle  if  the 
coefficient  is  negative).  For  example,  in  the  interaction 
diagram  above,  the  species  x  has  a  net  coefficient 
(f-l)/2  for  the  stabilizer,  because  f  >,  1.  This  is  valid 
because  the  matrix  elements  with  the  same  parameters 
cannot  be  distinguished  when  the  determinant  is  expanded. 

The  interaction  diagram  for  the  extreme  current  iM  becomes 


■- 
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(2.60) 


and  the  diagonal  elements  of  M  are:  -  -)  x ,  -fy ,  -2z. 

The  2  and  3-cycles  in  this  interaction  diagram 
correspond  to  the  product  of  matrix  elements  ^21 

and  Mj-2  M23  M3l  ^rom  respectively .  This  can  be 

deduced  by  considering  the  arrow 


i  . 


The  matrix  element  is  then  M. .  ,  where  Mn.  =  ax,  .  =  by 

ID  ll  2i  J 


ID 


2  j 


and  M_ .  =  cz  and  a,  b,  c  are  some  coefficients  which  are 
3d 

a  function  of  the  stoichiometric  parameter  f  in  some 
cases.  Note  that  the  components  of  £  are  assumed  to  be 
defined  in  the  order:  £  =  (X-X  ,  Y-Y  ,  Z-Z  )  as  in 

rU  O  O  O 

section  2.2.  In  conclusion,  the  cycles  which  are  found 


in  the  interaction  diagrams  of  a  network  give  directly 


. 


the  responsible  product  of  matrix  elements  in  the  inter¬ 
action  matrix. 

From  this  discussion,  it  is  clear  that  the  inter¬ 
action  diagram  above  is  identical  to  the  matrix  M  .  In 
general,  the  interaction  matrix  for  any  network  can  be 
expressed  as  a  sum  of  the  interaction  diagrams.  Since 
the  interaction  diagrams  correspond  to  the  extreme 
currents,  it  is  possible  to  derive  the  interaction  matrix 
by  formally  reducing  the  extreme  current  diagrams  to  the 

corresponding  interaction  diagrams.  This  is  done  by 

31 

applying  a  set  of  rules  given  by  Clarke: 

1)  A  solid  arrow  is  drawn  to  i  from  j  whenever  j  is 
a  reactant  and  i  is  a  product  of  the  same  reaction. 
The  number  of  barbs  on  the  new  arrow  equals  the 
number  of  left  feathers,  f  ,  times  the  number  of 
barbs  on  the  reaction  arrow. 

2)  A  dashed  arrow  is  drawn  to  i  from  j  whenever  j 
and  i  are  reactants  in  the  same  reaction.  The  number 
of  barbs  equals  the  product  of  f  at  j  and  the  total 
number  of  feathers  at  i. 

3)  The  value  of  the  stabilizer  of  each  reactant  is 
the  sum  over  all  outgoing  reactions  of  the  product 
of  the  numbers  of  total  feathers  and  f^. 

In  the  case  of  variable  stoichiometric  parameters, 
as  in  the  present  analysis,  the  arrows  of  the  interaction 


. 

. 
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diagrams  arc  labelled  instead  of  using  barbs,  for  con¬ 
venience.  Applying  Clarke's  rules  to  the  extreme  current 
diagrams  is  exactly  the  same  as  applying  the  procedure 

used  in  section  2.2  to  find  M,  from  the  extreme  current 

=A 

and  gives  an  alternate  method  of  deriving  the  inter¬ 
action  matrix  from  the  extreme  currents. 

Several  theorems  have  been  given  by  Clarke  that 
distinguish  the  types  of  cycle  diagrams  which  can 
destabilize  chemical  networks.  The  ones  pertinent  to 
this  analysis  are  the  following: 

1)  Cycles  with  even  numbers  of  dashed  arrows  (inhibitors) 
make  negative  contributions  to  ou.  (Theorem  4, 

Reference  19) . 

The  k-cycles  with  even  numbers  of  inhibitors  are  called 
; positive  feedback  k-cycles  or  PF  k-cycles ,  ’  If  a  positive 
feedback  k-cycle  contains  inhibitors  then  the  abbreviation 
PFI  k-cycle  will  be  used.  For  example,  in  the  inter¬ 
action  diagram  for  i^  given  in  (2.60),  a  PFI  2-cycle 
can  be  seen. 

2)  Cycles  with  an  odd  number  of  inhibitors  are  called 
negative  feedback  k-cy cles ,  or  NF  k-cycles.  The 
conditions  under  which  a  NF  3-cycle  can  destabilize  a 
network  is  given  in  the  following  theorem: 

A  NF  3-cycle,  which  passes  through  a  reactant  whose 
stabilizer  in  the  same  extreme  current  is  exactly  cancelled 


1 


\  3 


by  autocatalysis  (PF  1-cycle) ,  destabilizes  the  network 

if  no  other  cycles  containing  only  parameters  occurring 

in  the  3-cycle  pass  through  the  autocatalytic  species 

except  one  2-cycle.  (Theorem  3,  Reference  11). 

A  1-cycle  which  exactly  cancels  a  stabilizer  at  a  species, 

is  called  a  critical  1-cycle.  Although  a  critical  1-cycle 

does  not  violate  the  conditions  for  qualitative  stability, 

in  the  presence  of  a  negative  feedback  3-cycle  instability 

can  result.  In  the  interaction  diagram  for  i  given 

above,  the  negative  feedback  3-cycle  which  is  present 

does  not  pass  through  a  critical  1-cycle  for  f  ^  1,  and 

hence  this  cycle  cannot  destabilize  in  this  case.  However, 

when  f  =  1,  a  critical  1-cycle  appears  at  x,  and  the 

conditions  of  the  theorem  above  are  satisfied.  In  the 

results  of  Chapter  3,  it  will  be  shown  that  a  negative 

vertex  appears  in  IIA  corresponding  to  the  negative  feed- 

back  3-cycle  in  i  ,  at  f  =  1 . 

Up  to  this  point  the  extreme  current  and  interaction 

diagrams  that  have  been  considered  correspond  to  a  single 

extreme  point  in  the  current  polytope  II  .  In  principle, 

c 

similar  diagrams  also  exist  for  every  other  point  in  II  . 

However,  since  every  non-extreme  point  in  II  can  be 

c 

expressed  as  a  non-negative  linear  combination  of  the 

extreme  currents,  only  the  extreme  current  diagrams  are 

needed.  Thus,  for  any  given  point  in  n  ,  determined  by 

c 

the  current  vector  j ,  there  is  a  sum  of  weighted  interaction 


■ 
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diagrams  which  determine  the  interaction  matrix  M. 

Since  each  element  of  M  is  given  by  a  weighted  sum 
of  diagrammatic  elements,  one  from  each  extreme  current, 
the  expansion  of  the  determinant  for  M  yields  the 
characteristic  polynomial  as  a  homogeneous  sum  of  diagrams. 
Each  stability  inequality  (2.40)  in  the  network  can  also 
be  written  as  a  homogeneous  sum  of  diagrams.  It  follows 
that  a  set  of  diagrams  corresponds  to  the  set  of  points 
in  each  exponent  polytope,  II.  Therefore,  the  identif ca¬ 
tion  of  the  destabilizing  terms  in  f (p) ,  corresponding 
to  negative  vertices  of  II,  can  be  made  through  the 
diagrams.  The  diagrams  allow  easy  recognition  of  unstable 
states  which  have  only  an  abstract  topological  relation¬ 
ship  in  the  network.  In  order  to  determine  the  conditions 
for  instability,  the  adjacent  positive  vertices  must 
also  be  found.  Unfortunately,  rules  have  not  been  found 

that  can  identify  all  the  cycles  associated  with  the 
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adjacent  positive  vertices,  and  computer  programs  are 
used  to  calculate  the  adjacent  vertices  at  present. 

However,  the  diagrams  are  still  useful  for  eliminating 
certain  extreme  currents  from  the  analysis,  such  as  those 
which  do  not  lie  adjacent  to  the  destabilizing  extreme 
currents  in  the  current  polytope.  The  diagrams  also  retain 
information  about  the  chemical  interactions  which  cause 
instability,  and  so  pinpoint  the  cause  of  chemical 
instabilities . 
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CHAPTER  THREE 


The  Stability  Analysis 
3.1.  Introduction 

In  this  chapter  the  results  of  the  stability  analysis 
performed  on  the  Oregonator  model  (1.6)  are  given.  The 
irreversible  model  is  considered  in  section  3.2  and  the 
reversible  model  in  section  3.3.  These  results  are  then 
compared  in  section  3.5.  The  results  of  the  analysis 
give  the  change  in  the  boundary  of  the  unstable  steady 
states  when  the  reactions  are  allowed  to  reverse  in  the 
network.  Several  generalizations  emerge  and  are  given  in 
section  3.6. 

The  coefficients  of  the  terms  in  the  stability 
polynomials  are  a  function  of  the  stoichiometric  parameter 
f.  The  corresponding  vertices  in  the  exponent  polytope 
will  also  have  coefficients  which  are  a  function  of  f. 

The  surface  structure  of  the  exponent  polytope  will  change 
only  for  values  of  f  at  which  a  coefficient  of  a  vertex 
vanishes.  This  causes  a  truncation  of  the  polytope,  and 
the  region  of  instability  may  be  altered.  The  changes 
in  the  boundary  of  the  unstable  states  with  f  is  discussed 
in  detail  in  sections  3.2  and  3.3.  It  is  interesting  to 
note  that  adding  reverse  reactions  does  not  alter  the 
important  singular  values  of  f  which  are  present  in  the 


76 


■ 


' 


\ 

' 


irreversible  model.  These  occur  at  f 


0.5  ,  1  and 


1  +  /27 

The  complete  set  of  extreme  currents  for  the  network 
defines  terms  in  the  stability  polynomials  f (p)  which 
are  called  pure  or  mixed  terms.  The  pure  terms  are 
derived  from  one  extreme  current  only,  while  the  mixed 
terms  arise  from  the  overlap  of  different  extreme  currents. 
If  any  of  these  terms  is  dominant ,  it  can  be  used  as  a  one 
term  approximation  to  g  (tt)  =  0  and  is  therefore  a  vertex 
of  the  corresponding  exponent  polytope.  A  dominant 
negative  term  is  called  a  destabilizing  term.  The  purpose 
of  the  analysis  is  to  look  for  a  connection  between  the 
destabilizing  terms  in  f (p)  and  the  resulting  boundary 
of  stability. 


. 


•• 


3.2. 


Irreversible  Model  Results 


The  theory  developed  in  the  previous  chapter  can 
be  applied  to  the  irreversible  model  (1.6)  without 
difficulty  due  to  the  small  number  of  parameters  involved 
The  results  of  this  analysis  give  the  background  for 
comparison  with  the  results  obtained  when  reversibility 
is  included  in  the  Oregonator.  In  fact,  the  results  of 
the  irreversible  model  are  contained  in  the  reversible 
model  analysis,  and  represent  the  special  case  when  the 
reverse  rate  constants  are  set  equal  to  zero. 

The  interaction  matrix  is  derived  in  the  manner  of 
section  2.2,  with  the  stoichiometric  parameter  f  included 
as  an  additional  variable.  The  singular  nature  of  the 
value  f  =  1,  necessitates  the  division  of  the  analysis 
into  two  parts:  f  ><:  1  and  f  £  1 .  There  are  only  two 
extreme  currents  in  each  region  of  f,  as  shown  in  figure 
3.1a.  In  the  region  f  ><  1  the  extreme  currents  B  and  C 
exist,  while  for  f  ^  1  the  extreme  currents  A  and  C  exist 
Note  that  at  f  =  1  current  A  equals  current  B  as  required 
by  the  continuity  of  the  solutions  at  f  =  1 .  The  extreme 
current  A  was  derived  explicitly  in  section  2.2,  and  its 
diagrammatic  representation  was  given  in  section  2.5. 

The  interaction  diagrams  corresponding  to  each 
extreme  current  are  given  in  figure  3.1b.  Note  that  all 
the  elements  of  the  interaction  diagram  for  the  extreme 
current  A  have  been  multiplied  by  two  to  avoid  fractions. 


. 
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Figure  3 . 1 

(a)  The  extreme  currents  for  the  irreversible  model. 

(b)  The  interaction  diagrams  for  the  irreversible  model. 


The  interaction  matrix  M  in  each  region  of  f  is  determined 
directly  from  the  sum  of  the  diagrammatic  elements  of 
the  interaction  diagrams  B  and  C  for  f  ^  1,  and  A  and 
C  for  f  >,  1 . 

When  the  Lienard-Chipart  stability  conditions  (2.40) 
are  applied  to  the  characteristic  polynomial  for  M,  the 
stability  inequalities  take  the  form 

>  0,  >  0,  ~  aia2  “  aoa3  >  ® 

where  =  1.  However  as  will  be  discussed  in  section 

3.2.2,  the  analysis  of  is  simplified  if  the  inequality 

a2  >  0  is  included.  The  destabilizing  terms  were  found 

to  occur  in  two  of  these  inequalities,  and  only. 

The  responsible  destabilizing  cycles  can  be  chosen 

directly  from  the  interaction  diagrams1  in  figure  3.1b, 

on  the  basis  of  the  theorems  given  at  the  end  of  section 

2.5  which  identify  the  possible  destabilizing  cycles  in 

a  network.  Since  only  two  extreme  currents  are  present 

in  each  region  of  f,  the  exponent  polytopes  are  readily 

plotted  for  the  irreversible  model,  and  the  conditions 

for  instability  can  be  written  directly  from  each  exponent 

polytope,  II  and  II  .  No  destabilizing  terms  occur  in 
a2  A2 

or  which  are  positive  for  all  values  of  f.  The 
destabilizing  terms  in  and  A2  will  now  be  discussed 
and  a  summary  given  in  section  3.2.3. 


' 
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3.2.1.  The  Cone  of  Negative 

The  stability  inequality  for  >  0,  is  a  polynomial 

which  is  homogeneous  of  order  2  in  the  powers  of  each  of 
the  parameters  h  and  j .  The  parameters  will  be  denoted 
here  as  j^,  j  ,  x,  y,  z  in  the  region  f  >  1,  and  as  j  , 
jc,  x,  y,  z  in  the  region  f  <  1,  where  x,  y,  and  z  are  the 
reciprocals  of  the  steady  state  concentrations  Xq,  Yq , 

ZQ,  respectively.  Each  term  in  the  polynomial  is  a 
sum  of  all  the  2-cycles  and  2-stabilizers  which  occur  in 
the  interaction  diagrams  with  the  same  parameter  values. 

The  value  of  each  diagrammatic  element  is  read  from  the 
interaction  diagram  in  figure  3.1b  according  to  the 
discussion  given  in  section  2.5. 

The  only  destabilizing  cycle  in  is  the  inhibitory 
positive  feedback  (PFI)  2-cycle  between  x  and  y,  which 
occurs  in  the  extreme  currents  A  and  B  as  shown  in  figure 
3.1b.  The  destabilizing  term  in  which  corresponds  to 
this  2-cycle  is  found  in  each  region  of  f  by  summing  the 
pure  stabilizer  and  PFI  2-cycle  diagrams 

•  •  _  _  _  _ 

v  ~  —  — 

present  in  each  of  the  extreme  currents  A  and  B.  The 
coefficients  of  this  destabilizing  term  in  each  region 


of  f  are  as  follows: 


■ 
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f  <  1  f  =  1  f  >  1 

jgxy  f(l-2f)  -1 

j^xy  -4  2 (f 2-2f-l) 

2  2  2 

Note,  current  A  was  multiplied  by  2,  so  that  j  =  2  j 

A  B 

at  f  =  1 .  These  terms  are  negative  in  the  region 

0 . 5  <  f  <  1  +  /2~T  At  f  =  0.5  and  f  =  1  +  /2~ the 

coefficient  of  the  destabilizing  term  vanishes.  The 

remaining  terms  in  are  positive,  and  so  >  0  for 

both  these  values  of  f.  Similarly,  when  f  <  0.5  or 

f  >  1  +  /27  no  negative  terms  occur  in  •  Therefore, 

there  is  no  unstable  region  in  II  for  f  ^  0.5  and 

a2 

f  »  1  +  /2T 

The  unstable  region  in  IIo^  can  be1  discussed  separately 

in  each  of  the  regions  0.5  <  f  ^  1  where  current  B  exists 

and  1  $  f  <  1  +  /2~~ where  current  A  exists.  In  the  region 

2 

0.5  <  f  <  1  the  only  destabilizing  term  is  j_.xy  which 

corresponds  to  a  vertex  of  the  exponent  polytope  IIc^* 

This  vertex  has  three  edges  with  adjacent  positive  vertices 

.2  2  .2 

corresponding  to  the  terms  j„xy,  jDxz  and  J^yz.  The  three 

C  B  B 

edges  determine  three  inequalities  which  are  the  approximate 
necessary  and  sufficient  conditions  for  the  negative 
term  to  be  larger  than  the  sum  of  positive  terms  on  each 
edge.  Hence,  the  asymptotic  boundary  of  the  unstable 
region  in  a ^  j-s  given  by  the  hyperplane 
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tt  *  (y9  -  y,  )  =  log  |  c+/c  | 

where  tt  =  (log  jA,  log  jQ/  log  j  ,  log  x,  log  y,  log  z), 

and  the  cone  of  negative  is  the  intersection  of  the 
associated  halfspaces 

TT  •  (y9  -  yn  )  >  log  I  cVc“  I 

as  defined  in  section  2.4.  In  this  notation  for  tt  ,  the 

/v 

2 

vertex  corresponding  to  the  term  j  xy  has  the  exponent 

J 3 

vector  y0  =  (0,  2,  0,  1,  1,  0)  .  tt  is  defined  in  this 

way,  although  =  0  for  all  f  <  1 ,  to  maintain  consistency 

with  the  labelling  of  the  reversible  model  currents  of 
section  3.3.  The  cone  of  negative  is  defined  for  the 
region  0.5  <  f  <  1,  by  the  instability  conditions  given 
in  Table  3.1a.  Note  that  the  first  instability  condition 
corresponds  to  an  edge  of  II  containing  the  midpoint 

0 i  ^ 

A 

(0,  1,  1,  1,  1,  0),  which  has  been  included  (for  further 

details  see  the  end  of  section  2.4). 

In  the  region  1  <  f  <  1  +  /2"7  the  only  destabilizing 
2 

term  is  j  xy,  which  has  the  exponent  vector  y9  =  (2,  0, 

A  A 

0,  1,  1,  0).  The  terms  corresponding  to  the  adjacent 

.2  .2  .2  .  . 

positive  vertices  are  jcxy,  j^xz  aRd  jAyz.  T^e  con<^lt:Lons 

for  instability  in  this  region  are  given  in  Table  3.1b. 

The  terms  in  a 2  in  the  regions  0.5  <  f  <  1  and 

1  <  f  <  1  +  /2— differ  only  by  the  current  component  or 

jg.  Substituting  one  of  these  currents  for  the  other  in 


-  - 


■ 
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Table  3 . 1 

The  instability  conditions  for  in  the  irreversible  model, 
where  tt  =  (log  j^,  log  log  log  x,  log  y,  log  z)  for 

f  in  the  range  (a)  0.5  <  f  <  1,  (b)  1  <  f  <  1  +  /2j  (c)  f  =  l. 


0.5  <  f  <  1  Vertex  (0 , 2 , 0 , 1 , 1 , 0)  f(l-2f) 


Adjacent  Vertex 

Coefficient 

Instability  Condition 

1. 

(0,0, 2, 1,1,0) 

f (l+2f) 

(2f-l)(jB/jc)2  >  (l+2f)  +  2(l+f)jB/jc 

2. 

(0,2, 0,1, 0,1) 

2(l-f) 

f (2£-l)y  >  2(l-f)z 

3. 

(0,2, 0,0, 1,1) 

2f 

(2f-l)x  >  2z 

(a) 

1 

<  f  <  1  +  /2 

Vertex  (2,0, 

0,1, 1,0)  2(f2-2f-l) 

Adjacent  Vertex 

Coefficient 

Instability  Condition 

1. 

(0,0, 2, 1,1,0) 

f (l+2f) 

1/ 2 | f 2— 2f — 1 [ ( j A/ j  c ) 2  >  f (1+2 f)  + 

2f(l+3f)jA/jc 

2. 

(2, 0,0, 1,0,1) 

4(f-l) 

| f 2=2f-l | y  >  2 (f-l)z 

3. 

(2, 0,0, 0,1,1) 

8f 

| f 2-2f-l | x  >  4f z 

(b) 

f 

=  1  Vertex 

(0,2, 0,1, 1,0) 

2 

Adjacent  Vertex 

Coefficient 

Instability  Condition 

1. 

(0,0, 2,1, 1,0) 

3 

>  3  + 

2. 

(0,1, 1,1, 0,1) 

6 

iBy  >  6jcz 

3. 

(0,2, 0,0, 1,1) 

2 

x  >  2z 

(c) 


■ 

' 


■ 


* 

. 


any  term  of  only  modifies  the  coefficient  for  these 

ranges  of  f.  As  a  result,  the  instability  conditions 

given  in  tables  3.1a  and  3.1b  are  the  same,  except  for 

the  coefficients  which  shift  the  asymptotic  hyperplane 

defining  the  cone  of  negative  a^. 

At  f  =  1  the  extreme  currents  A  and  B  are  equal . 

It  follows  that  the  instability  conditions  for  each  region 

of  f  are  equal,  apart  from  the  factor  of  two  which 

multiplies  current  A.  However,  a  critical  1-cycle 

occurs  at  x  in  the  currents  A  and  B  when  f  =  1,  and  as 

2  2 

a  result  the  term  j^xz  (=  j^xz)  vanishes,  since  no  2-cycle 

passes  through  x  and  z.  This  truncates  the  exponent 

polytope,  allowing  a  different  exponent  vector  to  become 

a  vertex  of  II  .  The  result  is  that  a  new  positive 

a2 

vertex  (0 , 1 , 1 , 1 , 0 , 1)  occurs  in  II  for,  f  =  1,  and  a  new 

a2 

edge  is  formed  with  the  negative  vertex  (0,2, 0,1, 1,0). 

The  instability  conditions  for  f  =  1  are  given  in  table 
3.1c. 

The  exponent  polytope  II  can  be  visualized  in 

a2 

(log  x,  log  y,  log  z)  space  by  plotting  the  exponent 
vectors  corresponding  to  constant  current  terms  in  the 
plane  perpendicular  to  the  vector  £  =  (1,1,1).  Such  a 
plot  has  been  made  in  figure  3.2  for  various  values  of 
f.  The  maximal  polytope  is  truncated  at  f  =  0.5,  1  and 
1  +  /2  as  expected. 
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(Oil) 


0 . 5  <  f  <  1 


f  =  1 


1  <  f  <  1  +  y/2 


Figure  3.2 

The  exponent  polytope  n.  for  the  irreversible  model  at 

a  2 

various  values  of  f.  The  vector  tt  =  (1,1,1)  is  normal  to 

/V 

the  face  of  constant  j . 
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3.2.2.  The  Cone  of  Negative 


The  cone  of  negative  is  a  union  of  the  cones 
arising  from  each  negative  vertex  and  its  edges  in  . 
Since  >  0,  the  product  is  negative  whenever 

0.2  is  negative.  Therefore,  the  diagonal  product, 
of  [±2  defines  a  polytope  whose  set  of  negative  cone: 

arise  from  the  division  of  C  by  the  positive  dominant 

a2 

31 

terms  of  ou  .  IK  .  has  a  maximum  volume  when  all 
1  diag 

possible  parameter  combinations  occur  in  the  product 

ala2*  T^e  vertices  °f  ^diag  are  caiiec^  maximum  overlap 

(MO)  vertices  which  are  associated  with  maximum  overlap 

diagrams  (MOD's),  and  correspond  to  the  dominant  terms 

in  the  product  (see  theorem  10,  reference  19).  If 

all  the  vertices  of  IK  have  MOD's  then  II.  has  the  same 

A2  .  A2 

surface  structure  as  In  this  case,  the  boundary 

of  the  unstable  region  of  II 2  is  given  by  the  instability 
conditions  defining  the  cone  C 

a2 

If  is  truncated  from  the  maximum  polytope 

because  MOD's  are  missing,  then  new  vertices  will  be 


formed  in  II 


These  vertices  may  arise  from  terms  in 


A 2  which  do  not  have  maximum  overlap,  such  as  -a^  terms. 

A  vertex  of  IK  which  does  not  have  maximum  overlap  is 

2 

referred  to  as  a  non  maximum  overlap  (NMO)  vertex  and 
the  associated  diagrams  are  called  non  maximum  overlap 
diagrams  (NMOD’s).  When  terms  in  A2  which  have  NMOD ' s 
are  vertices  of  II  ^  the  cone  of  negative  A 2  may  destabiliz 


■ 
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regions  of  parameter  space  not  included  in  C  .  It 

a  2 

follows  that  the  only  vertices  that  need  to  be  considered 

in  II.  are  the  negative  NMO  vertices  which  occur  when 
2 

II.  is  truncated  from  its  maximum  configuration. 

A2 

The  presence  of  a  negative  NMO  vertex  in  IT,  is 
guaranteed  by  theorem  2  in  section  2.5,  which  identifies 
the  necessary  destabilizing  cycle.  According  to  this 
rule,  a  negative  feedback  (NF)  3-cycle  which  passes 
through  a  critical  1-cycle  is  a  possible  source  of  instab¬ 
ility  in  The  theorem  holds  because  the  NF  3-cycle 

acquires  a  minus  sign  in  the  expansion  of  the  determinant 
£±21  and  the  critical  1-cycle  eliminates  parameter  com¬ 
binations  in  a^a2  that  could  possibly  cancel  the  3-cycle. 

By  inspection  of  the  interaction  diagrams  in  figure  3.1b, 

i 

a  NF  3-cycle  is  found  in  currents  A  and  B.  However,  a 

critical  1-cycle  occurs  only  at  f  =  1 ,  when  the  coefficient 

of  the  x-stabilizer  vanishes  in  both  currents  A  and  B. 

The  destabilizing  term  in  A2  at  f  =  1,  which  corresponds 

.  3 

to  the  negative  NMO  vertex  found  m  n.  ,  is  ;j_.xyz.  The 

A  2  B 

diagrams  which  contribute  to  this  term  from  the  product 
ala2  w^en  f  differs  from  one  are  shown  in  figure  3.3.  The 
coefficient  of  this  NMO  term  is: 

j^xyz  2f  (2-3f )  f  <  1 

a 

j^xyz  16f  ( f — 2 )  f  >  1 

-2  f  =  1 


/ 


\ 
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Figure  3.3 


The  diagrams  contributing  to  the  NMO  vertex  in  II 


90. 


Note  that  at  f  =  1 .  The  coefficient  of  this 

term  is  negative  throughout  the  region  2/3  <  f  <  2. 

However,  within  this  range  of  f  the  polytope  is 

truncated  from  its  maximum  configuration  only  at  f  =1 , 

allowing  the  NMO  term  to  possibly  destabilize  additional 

regions  of  II.  . 

z 

The  truncation  of  the  poly  tope  II.  at  f  =  1  occurs 

2 

because  the  critical  1-cycle  at  x  in  current  B  (=A) 

3  2  3  2 

causes  the  maximum  overlap  (MO)  terms  jDx  y,  j^x  z  and 

r$  Jd 

.32. 

jbxz  ,  which  correspond  to  vertices  of  II ^  ,  to  vanish. 

Each  of  these  terms  can  be  identified  with  the  product 

ala2*  Since  >  0  for  f  /  1  (the  x-stabilizer  does  not 

vanish),  II  , .  will  be  maximal  when  II  is  maximal.  From 

diag  a2 

the  discussion  in  section  3.2.1,  II'  is  only  truncated 

2 

when  f  =  0.5,  1  and  1  +  /2T  The  case  f  =  1  has  already 

been  discussed,  where  the  truncation  in  IT.  is  due  to 

.2  2 

the  a2  term  jBxz  vanishing  as  a  result  of  the  critical 

1-cycle  at  x.  When  f  =  0.5  and  1  +  /tT,  the  coefficients 

2  2  . 

of  the  terms  jBxy  and  3^xy  vanish  respectively.  This 


•  3  2 

causes  both  of  the  corresponding  MO  terms  j  x  y  and 
3  2 

j  xy  to  vanish  in  A2,  truncating  the  polytope 
Since  there  is  no  critical  1-cycle  passing  through  the 
NF  3-cycle  for  these  values  of  f,  no  negative  NMO  vertex 
exists  to  cause  possible  additional  destabilization  in 


.3 


.3 


II  ^  .  In  fact  the  jDxyz,  j^xyz  terms  are  positive  and 


B 


A 


lie  on  an  edge  of  II 


A, 


for  f 


0 . 5  and  f 


1  +  /2 


\ 


is  sketched  for 


respectively.  The  exponent  polytope  II 


A. 


several  values  of  f  in  figure  3.4.  Once  again  the 
vector  (1,1,1)  in  (log  x,  log  y,  log  z)  space  is  perpen¬ 
dicular  to  the  face  of  constant  j . 

The  instability  conditions  for  the  one  NMO  negative 
vertex  at  f  =  1  are  given  in  Table  3.2  where  tt  is  the 
same  as  previously  defined  in  section  3.2.1.  The  instab¬ 
ility  conditions  define  hyperplanes  with  the  positive 
adjacent  vertices  which  are  identical  to  those  defined 

by  C.  for  f  =  1.  Therefore,  C  has  not  been  extended 
a2  a2 
by  the  cone  of  negative  due  to  the  negative  NMO  vertex 

present  at  f  =  1  in  this  model. 


' 


■ 
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Figure  3.4 


The  exponent  polytope  II 


A, 


for  the  irreversible  model. 


shown 


for  various  values  of  f 
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Table  3 . 2 

The  instability  conditions  for  the  negative  NMO  vertex  in 
A2  at  f  =  1,  where  =  (log  j  ,  log  jB,  log  j  ,  log  x,  log  y, 
log  z ) . 


NMO  Vertex 

(0,3, 0,1, 1,1, 

)  -2 

Adjacent  Vertex 

Coefficient 

Instability 

Conditions 

1. 

(0,2,1, 1,0,2) 

12 

V  > 

6jcz 

2. 

(0,3, 0,0, 1,2) 

4 

x  > 

2z 

3. 

(0,3,0,1,2,0) 

-1 

2z  > 

y 

4. 

(0,2, 1,2, 1,0) 

-3 

2V  > 

3jcx 

3.2.3.  Summary 


In  the  irreversible  model,  one  negative  term  occured 
in  a  stability  inequality  according  to  the  following  table 

0 . 5<f  <1  f=l  1< f  <  1+/2  f  >l+/2 

-  + 

+  +  +  + 

+  -  +  + 

The  responsible  destabilizing  cycle  and  the  extreme 
current  in  which  the  cycle  originated,  for  each  destabil¬ 
izing  term  in  the  table  above,  is  summarized  in  Table  3.3. 
One  pure  destabilizing  term  occurs  in  and  destabilizes 
the  network  for  0.5  <  f  <  1  +  /2T  At  f  =  1,  one  pure 
NMO  destabilizing  term  occurs  in  due  to  the  negative 

feedback  3-cycle,  but  does  not  extend  the  boundary  of 
the  unstable  states  defined  by  Cc^. 

The  asymptotic  boundary  of  the  unstable  steady  states 
in  the  network  is  given  by  the  three  instability  condi¬ 
tions  which  are  listed  in  Table  3.1  for  each  region  of 
f.  For  f  j*  1,  two  of  the  instability  conditions  define 
hyperplanes  in  the  log  of  reciprocal  concentration  space 
independently  of  the  currents.  These  hyperplanes  have 
been  plotted  in  Figure  3.5  for  the  case  f  =  1.1.  If  the 
value  of  f  is  changed,  then  the  planes  defining  the 
asymptotic  boundary  are  shifted  parallel  to  the  axis 


0<f$0. 5 


a. 


a. 


+ 

+ 

+ 


/ 
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Table  3.3 


Summary  of  the  cycles  causing  instability  in  the  irreversible 
model. 


Destabilizing  cycles 

The  extreme 

Unstable  range 

currents  involved 

of  f 

a 

pure  A 

i  f  <  i  +  rr 

pure  B 

0.5  <  f  <:  1 

none 

N.  * 

• 

a2(nmo)  / 

pure  A 

X 

f  =  1 

pure  B 

f  =  1 

* 


A  critical  1-cycle  is  necessary  in  addition. 
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Plot  of  the  irreversible  model  hyperplanes  at  f=l.l 
unstable  region  is  designated  by  U. 


The 
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intercepted  in  the  positive  orthant.  The  result  is  that 

the  "wedge"  of  unstable  parameter  space  is  also  shifted. 

The  relationship  between  these  asymptotic  hyperplanes 

and  the  exponent  polytopes  II  and  IL  is  indicated  by 

a2  2 

the  triangle  and  hexagon  respectively,  outlined  perpen¬ 
dicular  to  the  vector  (1,1,1)  in  Figure  3.5.  The  exponent 

polytopes  II  and  nA  of  Figures  3.2  and  3.4  are  reduced 

2  2 

to  two  dimensions  by  taking  the  current  dependency  of 

each  vertex  into  the  coefficient  of  that  vertex,  and 

then  summing  these  coefficients  over  the  exponent  vectors 

with  the  same  reciprocal  concentration  parameters  (i.e., 

the  points  of  II  or  IIA  which  lie  on  a  line  parallel 

2  2 

to  the  vector  (1,1,1)  are  represented  by  one  point  of 

the  triangle  or  hexagon  with  a  coefficient  which  depends 

on  the  currents  and  f ) .  Figure  3.5  is  interesting 

because  it  shows  how  the  hyperplanes  defined  by  CL  can 

2 

become  extended  into  CA  .  From  Figure  3.5,  it  is  not 

2 

surprising  that  CA  did  not  extend  the  boundary  defined 

2 

by  C  .  The  dashed  lines  on  the  (log  x,  log  z)  plane 
2 

give  the  position  where  the  asymptotic  planes  of  the 
unstable  region  in  the  f  =  1.1  case  cut  the  (log  x,  log  z) 
plane.  The  projection  of  the  hexagon  onto  the  (log  x, 
log  z)  plane  is  also  outlined  in  Figure  3.5,  and  can  be 
compared  to  Figure  2.1  which  corresponds  to  the  case 


jc  =  °'  3a  =  1  and  y  = 


1. 


’ 

If  f;!  :1  1  •'  .  1 


In  section  2.2  the  steady  state  represented  by  the 
overall  reaction, 

°f:  1+f  ~A  +  ITf  ~C  '  f  >  1  ' 

was  given.  To  test  the  stability  of  the  steady  state 
represented  by  this  combination  of  the  extreme  currents, 
the  ratio  of  jA/jc  which  must  be  violated  for  instability 
is  needed.  The  criterion  for  an  unstable  ratio  of  these 
currents  is  determined  by  the  instability  condition 

1/2 | f 2-2f-l [ ( jA/jc) 2  >  f(l+2f)  +  f (l+3f) (jA/jc)  (3. 

given  in  Table  3.1b  for  (f  >  1).  The  inequality  (3.1) 
implies  that  j A/ j c  must  be  greater  than  the  positive 
root  of  the  quadratic  equation  obtained  by  substituting 
an  equal  sign  for  the  inequality.  For  the  case  f  =  1.1 
this  equation  yields  jA/j^  >  5.4,  while  the  linear  combin¬ 
ation  0f  above  has  j A/ j c  =  20  for  f  =  1.1.  Therefore, 

Of  is  a  possibly  unstable  steady  state  for  f  =  1.1.  For 
any  particular  steady  state  concentration  of  the  species, 
the  remaining  instability  conditions  in  Table  3.1b  must 
be  tested,  before  instability  can  be  determined.  In  the 
case  f  =  1.5  the  instability  condition  (3.1)  gives 
3A/jc  >  10  *1  while  °f  yields  j  /j  =  4.  This  means  that 
Of  represents  a  stable  steady  state  for  f  =  1.5.  It 
follows  that  Of  is  possibly  unstable  only  for  values  of 

In  addition,  instability  condition  (3.1) 


f  close  to  one. 


\ 


' 
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implies  that  as  f  increases,  the  current  necessary  in  the 
extreme  current  A(f  >  1)  for  instability  also  increases. 

The  instability  conditions  give  the  requirements 
on  the  parameters  for  unstable  steady  states.  As  such 
they  do  not  guarantee  instability  will  occur,  only  that 
the  parameters  can  be  adjusted  so  that  the  steady  state 
will  become  unstable.  The  actual  rate  constants  may 
not  allow  a  portion  of  the  unstable  region  to  be  attained 
in  a  network. 

The  currents  of  the  irreversible  model  can  be  written 
in  terms  of  the  rate  constants  by  using  equations  (2.33) 
and  the  extreme  currents: 


i  =  (1+f>  1  0  1) 

-A  K  *5  '  o  ,  1/  U,  i; 


2  '  2 


i_  =  (f,  0,  1,  ^  ,  1) 

^  L.  Z 


which  are  valid  in  the  region  f  >  1.  The  following 
relations  are  then  obtained: 


kl  =  K1A  =  1/Yo  3a  +  f jc} 


(1+f)  . 

o  o  ‘  2  JA 


k2  =  K2  =  1/X„Y„  j*} 


=  K^B  =  1/XQ  {jft  +  j„} 


k4  -  K4  “  1/Xo  {i^  V 


k5  “  Ks  =  1/ZQ  {jA  +  jc} 


' 
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Using  the  second  and  fourth  relations  above,  the  instability 
condition  (3.1)  becomes 

V*c  =  K2yo/K4Xo  >  5'4 

at  f  =  1.1.  Since  the  steady  state  concentrations  are 
functions  of  the  rate  constants  also,  the  relations 

Z  =  K-.BX  /Kc 
o  3  o  5 

Yq  =  fK3BXo/(K1A  +  K2Xq)  (3.2) 

-2K2K4X^  +  XQ ( (1-f )K2K3B-2K1K4A)  +  (l+f)K1K3AB  =  0, 

which  are  obtained  from  the  solution  of  the  kinetic 
equations  at  steady  state,  can  be  used  to  eliminate  the 
concentration  parameter  from  the  •  instability  conditions. 

At  f  =  1  the  instability  condition  x  >  2z  in  Table 
3.1c  can  be  written  as 

K3BXq/K5  >  2Xq  (3.3) 

where  the  equation  for  ZQ  in  (3.2)  has  been  used  in 

addition  to  the  definitions  x  =  1/X  ,  y  =  1/Y  .  This 

o  o 

equation  can  be  used  to  give  the  upper  limit  on  for 

which  instability  occurs  at  f  =  1 .  In  order  to  compare 

24 

this  value  to  the  one  obtained  by  Field  and  Noyes, 

-2  7  -1  -l 

the  values  B  =  6  x  10  M  and  K3  =  8  x  10  M  sec  used 
in  their  analysis  are  substituted  into  equation  (3.3). 

The  factor  of  2  in  front  of  Z  was  absent  in  step  (5)  of 


-  -<r,  .  t  ' 


their  model  (compare  with  (1.6)),  and  so  the  factor  of 
2  in  equation  (3.3)  is  also  absent.  Finally,  this  yields 

K5  <  K3B  =  480  sec"1 


for  unstable  states  in  the  Oregonator  model  at  f 
The  value  obtained  by  Field  and  Noyes  was  = 


1. 

-1 

sec 
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3.3.  The  Reversible  Model 


In  this  section  the  boundary  of  the  unstable  steady 
states  for  the  reversible  model  is  given  in  a  parallel 
fashion  to  the  preceding  irreversible  case.  The  kinetics 
used  for  each  step  of  the  reversible  model  is  determined 
by  the  following  table,  where  the  indicated  kinetic  order 
is  with  respect  to  each  species:^ 


forward  kinetic 

reverse  kinetic 

order 

order 

1.  A  + 

Y 

X  +  P 

1st 

1st 

2.  X  + 

Y 

2P 

1st 

1st 

3.  B  + 

X 

2X  +  2Z 

1st 

2nd  in  X,  1st  in 

4.  2X 

Q 

2nd 

1st 

5.  2Z 

"V. 

- -  fY 

1st , 

1st 

These  reactions  yield  thirteen  extreme  currents  in 
each  region  of  f,  which  are  given  in  figure  3.6.  The 
extreme  currents  labelled  C,  C- ,  D,  D_ ,  1,  2,  3,  4  and  5 
are  valid  for  all  values  of  the  stoichiometric  parameter 
f.  The  four  remaining  extreme  currents  A,  A  ,  B  and  B 
differ  in  the  regions  f  <  1  and  f  >  1,  but  are  equivalent 
at  f  =  1.  Note  A  =  B  and  A  =  B  at  f  =  1. 

In  a  reversible  network,  the  extreme  currents  can 
be  divided  into  two  groups:  the  equilibrium  extreme 
currents  (EEC's)  and  the  non-equilibrium  extreme  currents 
(NEEC’s).  The  EEC's  consist  of  a  reaction  and  its  reverse 


", 


■ 
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f  <  1 


f  >  1 


all  f 


Figure  3.6 


The  extreme  currents  of  the  reversible  model. 
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only,  and  are  identified  with  the  extreme  currents  labelled 
1  to  5  (in  figure  3.6).*  The  NEEC ' s  contain  two  or  more 
reactions,  and  no  two  reactions  in  a  NEEC  are  the  reverse 
of  each  other.  The  NEEC ' s  in  figure  3.6  are  labelled  A 
to  D. 

The  EEC's  represent  the  steady  state  current  in  a 

reaction  when  the  forward  and  reverse  rates  are  equal, 

or  in  other  words  the  EEC's  satisfy  detailed  balance. 

4  6 

Shear  has  proved  that  chemical  systems  which  are 
decomposible  into  elementary  reactions  satisfying  mass 
action  kinetics  always  have  stable  steady  states  when 
detailed  balance  holds.  From  this  it  follows  that  any 
positive  linear  combination  of  the  EEC's  which  have  mass 
action  kinetics  must  be  stable  for  all  values  of  the 
parameters.  Only  the  EEC's  1,  2  and  4  have  mass  action 
kinetics  in  the  reversible  model  for  f  ^  1.  However, 
the  remaining  EEC's  3  and  5  do  not  have  any  distinguish¬ 
ing  behavior  as  a  result  of  their  non-mass  action  kinetics, 
and  so  need  not  be  discussed  separately  in  this  model. 

The  thermodynamic  equilibrium  steady  state  will  be  some 
positive  linear  combination  of  the  EEC's  in  the  network. 

Such  a  steady  state  will  rule  out  any  net  flow  of  matter 

* 

Because  of  difficulties  which  may  arise  in  non-mass  action 
systems,  future  articles  will  use  a  stronger  definition  of 
EEC  which  restricts  the  reactions  that  make  up  an  EEC  to 
those  having  mass  action  kinetics. 


- 

- 


' 
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in  a  closed  reaction  sequence  (NEEC) ,  since  the  flow  in 
the  forward  direction  is  equal  to  the  flow  in  the  reverse 
direction. 

Far  from  equilibrium,  a  steady  state  can  exist 
without  the  reactions  satisfying  detailed  balance, 
although  a  flow  of  matter  through  the  system  is  necessary 
to  maintain  the  steady  state  (open  system) .  A  NEEC 
represents  a  possible  steady  state  in  which  the  reverse 
of  the  reactions  present  in  the  NEEC  have  zero  rate  con¬ 
stants,  and  thus  represents  a  possible  steady  state  far 
from  equilibrium.  In  the  special  case  where  all  the 
reverse  rate  constants  are  set  equal  to  zero,  the  NEEC ' s 
which  exist  are  those  found  in  the  corresponding  irrevers¬ 
ible  network. 

The  current  polytope  II  for  the  reversible  model 

o 

lies  in  a  six  dimensional  subspace  of  the  seven  dimensional 

space  in  which  v°  lies,  as  the  application  of  the  theory 

in  section  2.2  shows.  Any  general  steady  state  of  the 

network  is  a  positive  linear  combination  of  the  extreme 

currents  and  is  represented  by  a  point  in  II  .  When  the 

steady  state  is  given  by  an  extreme  current,  the  system 

is  represented  by  an  extreme  point  of  II  ,  and  the  reactions 

c* 

not  included  in  the  extreme  current  have  zero  velocity. 

Two  extreme  points  are  adjacent  in  II  when  they  have  the 

o 

maximum  number  of  reactions  in  common,  since  this  is 
equivalent  to  having  the  maximum  number  of  components  of 


•• 


g 
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the  corresponding  basis  vectors  (extreme  currents)  in 

common.  The  NEEC ' s  A  and  B  have  three  out  of  four  possible 

reactions  in  common  as  do  the  NEEC's  B  and  C,  and  so 

these  NEEC 1  s  lie  on  a  common  face  of  II  .  Similarly  the 

c 

NEEC's  A-,  B~  and  C  lie  on  a  common  face  of  II  .  However 
none  of  the  extreme  currents  A  ,  B  and  C  are  adjacent 
to  A,  B  or  C.  The  NEEC  D  has  two  out  of  three  possible 
reactions  in  common  with  the  extreme  currents  C  and  A 
and  so  is  adjacent  to  these  NEEC ' s ,  while  the  NEEC  D  is 
adjacent  to  C~*  and  A  when  f  <  1.  In  the  region  f  >  1, 
the  NEEC  D  is  adjacent  to  C  and  B  ,  while  D  is  adjacent 
to  C“  and  B.  The  EEC's  represent  five  extreme  points  in 
a  four  dimensional  subspace  of  II  which  is  a  simplex 
since  each  EEC  is  adjacent  to  every  other  EEC.  The 
relationship  of  the  extreme  currents  to  each  other  in  II 

c 

is  shown  in  figure  3.7. 

On  the  basis  of  the  destabilizing  cycles  given  in 

the  theorems  of  section  2.5,  the  network  is  unstable 

when  the  currents  in  the  network  are  close  to  the  NEEC's 

A  or  B  (the  AB  edge  of  II  )  .  The  extreme  currents  C,  D  , 

c 

1,  2f  3,  4,  and  5  are  adjacent  to  either  A  or  B  in  II  , 
while  the  remaining  extreme  currents  A  ,  B  ,  C  and  D  are 
not  adjacent  to  either  A  or  B.  Therefore,  the  extreme 
currents  A-,  B~ ,  c”  and  D  were  eliminated  from  this 
analysis.  The  effect  of  the  deleted  NEEC's  on  the  net¬ 
work  will  be  discussed  in  section  3.4.  Even  with  this 
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Figure  3 , 7 

The  relationship  between  the  extreme  currents  of  the 
reversible  model  in  the  current  polytope  IIc. 
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reduction  in  the  number  of  extreme  currents,  the  analysis 
still  requires  thirteen  parameters  including  the  reciprocal 
steady  state  concentrations  x,  y,  z  and  the  stoichiometric 
parameter  f.  The  set  of  extreme  currents  chosen  form  a 
subset  of  the  complete  set,  and  will  serve  to  illustrate 
the  effect  of  adding  reverse  reactions  to  an  unstable 
irreversible  network. 

The  interaction  diagrams  corresponding  to  the  extreme 
currents  A,  B,  C,  D  ,  1,  2,  3,  4  and  5  are  given  in 
figure  3.8  for  both  regions  of  f.  The  NEEC  A  has  been 
multiplied  by  two  to  remove  fractions  as  in  the  irrevers¬ 
ible  model.  Once  again  the  Lienard-Chipart  conditions 

>  0,  >  0,  A2  >  0 

are  used  in  addition  to  the  inequality  >  0.  At  least 
one  negative  term  occurred  in  these  inequalities  accord- 

f  <  1  f  =  1  f  >  1 


As  discussed  in  the  previous  section,  only  the  NMO 
destabilizing  terms  are  considered  in  since  only 

these  terms  can  cause  appreciable  extension  of  the 
unstable  region  as  defined  by  the  union  of  the  negative 


ing  to  the  following  table: 


A2  (NMO) 


■ 


F<  1 


F>1 


2  y-==:::^x 
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/  /  1 


2z 
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Figure  3.8 


The  interaction  diagrams  used  in  the  reversible  model 
analysis . 


A  quick  comparison  to  the  table 


cones  of  and  a^. 
given  in  section  3.2.2  for  the  irreversible  case  shows 
that  the  negative  terms  in  (f  >  1)  and  A2  (f  <  1)  are 
new.  The  destabilizing  terms  found  in  each  region  above 
will  now  be  given,  together  with  the  instability  condi¬ 
tions  derived  from  the  edge  approximation.  A  summary  of 
the  results  appears  in  section  3.3.4. 


! 
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3.3.1.  The  Cone  of  Negative 

The  possibly  destabilizing  terms  in  can  be  selected 

from  the  interaction  diagrams  in  figure  3.8  by  summing 

the  inhibitory  positive  feedback  (PFI)  2-cycles  and 

stabilizer  diagrams  as  was  done  for  the  irreversible 

model  shown  in  section  3.2.1.  However,  the  large  number 

of  current  parameters  in  the  reversible  model  necessitated 

the  use  of  computer  programs,  developed  for  finding  the 

vertices  and  edges  of  exponent  polytopes  in  higher 

4  5 

dimensional  spaces,  before  the  instability  conditions 
were  determined.  In  addition,  programs  which  calculate 
the  terms  in  the  stability  inequalities  were  used  to 
expand  the  determinant  in  the  characteristic  equation 
for  y.42 

i 

Four  destabilizing  terms  exist  in  a.^,  two  of  which 
are  pure  terms  derived  from  the  NEEC's  A  and  B.  The 
remaining  two  are  mixed  terms  involving  either  current  A 
or  B.  The  destabilizing  cycle  in  each  term  is  the  PFI 
2-cycle  which  occurs  in  both  NEEC  A  and  B.  Table  3.4 
tabulates  the  coefficients  and  the  unstable  range  of  f 
for  each  destabilizing  term  in  c^/  where  the  correspond¬ 
ing  exponent  vectors  are  defined  by 

*  =  dog  jA,  log  jB,  log  jQf  log  jD~,  log  log  j2# 

log  j3,  log  j4,  log  j5,  log  x,  log  y,  log  z) . 

The  fifth  exponent  vector  in  Table  3.4  is  a  vertex  of 

the  exponent  polytope  n  only  at  f  =  0.5  and  1  4-  /2~7 

a2 
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2  2 

and  is  a  midpoint  of  the  j  -  j  edge  for  other  values 

B 

of  f.  At  f  =  0.5  the  pure  destabilization  j^xy  vanishes, 

B 

while  at  f  =  1  +  /2  the  pure  destabilizing  term  jfxy 

vanishes,  truncating  the  polytope  II  and  allowing  the 

a2 

midpoint  corresponding  to  the  term  j,jDxy  to  become  a 

A  B 

vertex  of  n  .  Note  from  Table  3.4  that  the  destabiliz- 
a2 

ing  terms  in  ol^  change  sign  at  either  f  =  0.5,  2  or  1  +  /2. 

The  a2  destabilizing  terms  all  contain  the  PFI 
2-cycle  which  passes  through  x  and  y  in  the  interaction 
diagrams  of  figure  3.8.  The  sum  of  this  2-cycle  and 
the  stabilizer  diagram  for  the  EEC  2  has  a  coefficient 
of  zero,  and  so  the  EEC  2  cannot  form  a  pure  destabiliz¬ 
ing  term  in  o^.  As  a  result,  the  mixed  destabilizing 

terms  j ^ j 2X^  anc^  ^b^2x^'  which  would  have  corresponded 

•  •  2  2  2  2 
to  midpoints  of  the  jA  -  j2  and  -  j2  edges  of  n 

respectively,  correspond  to  vertices  of  II  instead. 

a2 

The  cone  of  negative  a 9,  C  ,  for  any  value  of  f  is 

z  a2 

the  union  of  the  cones  C  resulting  from  each  negative 

a2 

vertex  in  II  .  Each  of  the  cones  C1  gives  the  region 
a2  a2 

of  dominance  of  the  corresponding  negative  vertex  i  in 

II  .  Since  the  coefficients  of  the  terms  in  a0  are,  in 
a2  ^ 
general,  quadratic  functions  of  f,  the  bounding  hyper¬ 
planes  of  the  cones  C1  will  shift  as  f  varies.  At  the 

a2 

values  of  f  where  a  coefficient  vanishes,  the  exponent 
polytope  becomes  truncated  and  new  conditions  for 
instability  appear,  corresponding  to  new  edges  of  II 

a2 


' 

' 


1 
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In  figures  3.9a  and  3.9b  the  adjacent  vertices  of 

the  pure  NEEC  negative  vertices  corresponding  to  the 
•  2  2 

terms  3^xy  and  3Bxy  are  shown  in  each  region  of  f.  The 

minimum  number  of  edges  occurs  in  the  region  f  >  1,  where 

the  exponent  polytope  is  maximal.  The  adjacent  vertices 

in  the  region  f  >  1  are  indicated  by  the  first  column 

in  figures  3.9a  and  3.9b.  When  a  vertex  vanishes,  the 

polytope  II  becomes  truncated  and  the  appearance  of  new 
a2 

vertices  occurs  in  a  regular  way.  The  adjacent  vertices 

in  figures  3.9a  and  b  which  vanish  at  either  f  ^  1  or 

f  =  1  are  indicated  by  a  star  and  a  double  star  respectively. 

The  new  vertices  which  then  appear  adjacent  to  the  pure 

negative  vertices  are  indicated  by  the  remaining  columns 

of  figures  3.9a  and  b.  The  interesting  feature  in  these 

figures  is  that  the  new  vertices  which  appear  are  all 

derived  from  one  of  the  starred  vertices  by  substituting 

one  of  the  current  components  of  the  starred  vertex  by 

the  current  component  of  each  of  the  other  extreme 

currents  in  the  network.  An  exception  occurs  when  the 

EEC  5  is  involved  in  two  vertices  which  both  vanish. 

In  the  region  f  >  1 ,  no  NEEC  has  a  critical  1-cycle 

and  the  exponent  polytope  II  is  maximal  unless  f  =  2, 

a2 

1  +  /2.  The  truncation  of  the  polytope  which  occurs  at 
f  =  2  and  f  =  1  +  /2  is  minimal  because  only  the  coefficient 
of  a  destabilizing  term  vanishes.  This  can  be  compared 
to  the  situation  in  the  irreversible  model  at  f  =  1  +  /2. 
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f  >  1 


f  $  1 


(200000000101)* - (110000000101)** 

(101000000101) 

(100100000101) 

(100010000101) 

(100001000101) 

(100000100101) 

(100000010101) 

[(100000001101)*]  - =*  (010000001101)** 

(001000001101) 
(000100001101) 
(000010001101) 
(000001001101) 
(000000101101) 
(000000011101) 


(100000001011) 


(100000100110)  same 

(100000010110)  same 

(100000001110)*  - ^  (010000001110)** 

(001000001110) 
(000100001110) 
(000010001110) 
(000001001110)  ' 
(000000101110) 
(000000011110) 


(200000000011) 

(020000000110) 

(002000000110) 

(000200000110) 

(100010000110) 

(100001000110) 


Figure  3.9a 

The  appearance  of  new  adjacent  vertices  as  f  varies  for 

2 

the  pure  destabilizing  term  j  xy,  which  corresponds  to 
the  vertex  (200000000110)  in  II 

“2 


* 


Terms  vanish  at  f  =  1  only. 


Terms  vanish  when  f  <  1. 


** 
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f  1  1 


f  =  1 


f 


1 


(020000000011)  same 

(020000000101) ** - ^ (011000000101) 

(010100000101) 

(010010000101) 

(010001000101) 

(010000100101) 

(010000010101) 

[ (010000001101) **] 


(200000000110) 
(002000000110) 
(000200000110) 
(010010000110) 
(010001000110) 
(010000100110) 
(010000010110) 
(010000001110)  * 


(000010001110) 

(000001001110) 

(000000101110) 

(000000011110) 


(001000001101) 

(000100001101) 

(000010001101) 

(000001001101) 

(000000101101) 

(000000011101) 


(010000001011) 


Figure  3 . 9b 

The  appearance  of  new  adjacent  vertices  of  the  pure 

2 

destabilizing  term  j  xy  corresponding  to  the  vertex 

a 

(020000000110)  in  II  .  The  ten  vertices  in  the  first 

a2 

column  are  adjacent  to  (020000000110)  for  f  <  1  and  f  >  1 

(except  special  cases) . 

** 

Terms  vanish  at  f  =  1  only  due  to  a  critical  1-cycle  at 


n  in  the  NEEC  B. 


. 
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When  f  >  1,  the  negative  vertices  in  II  have  only  ten 

a2 

adjacent  vertices. 

In  the  region  f  <  1,  the  NEEC  A  has  a  critical 

1-cycle  at  x.  The  resulting  truncation  of  the  exponent 

polytope  is  severe  because  no  2-cycle  exists  between  x 

and  z  in  the  NEEC  A,  and  the  vertex  corresponding  to  the 

.2 

term  j^xz  vanishes  causing  possibly  eight  new  adjacent 

vertices  to  form  with  the  vertex  corresponding  to  the 
.  2 

term  j  xy.  However,  one  of  the  new  adjacent  vertices 

involves  the  EEC  5  which  does  not  have  an  x-stabilizer 

for  any  value  of  f.  As  a  result,  this  vertex,  represented 

by  a  square  bracket  in  figure  3.9a,  vanishes.  In  addition, 

one  other  adjacent  vertex  present  at  f  >  1  involving 

the  EEC  5  vanishes  for  f  <  1.  Hence  a  further  truncation 

of  II  occurs  and  fifteen  new  edges  appear  as  shown  in 
a2 

figure  3.9a. 

When  f  =  1,  the  NEEC  B  also  has  a  critical  1-cycle 

at  x.  As  a  result,  three  of-  the  adjacent  vertices 

present  at  f  <  1  vanish  in  figure  3.9a,  but  no  new  adjacent 

vertices  form.  This  may  be  attributed  to  the  severe 

truncation  of  II  already  present  at  f  <  1  about  the 

a2 

...  .2 

vertex  corresponding  to  the  pure  destabilizing  term  3  xy. 

A  further  truncation  of  the  polytope  at  f  =  1  involves 

vertices  with  current  components  from  the  NEEC  B.  The 

.  .  .  .  2 

vertex  corresponding  to  the  pure  destabilizing  term  j  xy 
has  nineteen  new  edges  which  replace  two  vanishing  edges 


' 

. 
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of  II  at  f  =  1  as  shown  in  figure  3.9b.  Actually  the 
2 

truncation  of  IT  at  f  =  1  is  artificial  since  the  NEEC '  s 

a2 

A  and  B  are  equal  at  f  =  1.  However,  the  truncation  of 

n  and  the  subsequent  appearance  of  new  positive  vertices 
2 

pictured  in  figures  3.9a  and  b  explains  the  changes  in 
the  instability  conditions  which  occur  for  the  pure  A 
and  B  destabilizing  terms  in  each  region  of  f. 

The  instability  conditions  for  each  destabilizing 
term  in  are  given  in  tabular  form  for  each  region  of 
f  in  which  they  apply.  Table  3.5  summarizes  which  of 
the  sets  of  instability  conditions  given  in  Tables  3.6 
to  3.10  are  applicable  in  each  region  of  f. 

An  examination  of  the  instability  conditions  for 
(*2  in  each  region  of  f  yields  several  important  results. 
First,  in  the  region  f  >  1  the  only  difference  in  the 
adjacent  positive  vertices  of  the  pure  A  and  B  negative 
vertices  is  that  current  A  replaces  current  B.  Similarly 
the  mixed  current  negative  vertices  have  adjacent 
vertices  which  differ  only  by  the  current  components  A 
and  B.  This  is  to  be  expected  since  the  interaction 
diagrams  for  the  NEEC ' s  A  and  B  in  the  region  f  >  1  differ 
only  by  the  coefficient  of  the  PFI  2-cycle. 

In  the  region  f  <  1  a  similar  parallel  amongst  the 
pure  and  mixed  negative  vertices  is  observed,  except  for 
the  additional  instability  conditions  which  appear  for 
the  pure  A  vertex. 
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Table  3.5 

A  guide  to  the  use  of  the  tables  of  instability  conditions 
for  a2  in  the  reversible  model:  the  region  of  f  in  which 
the  instability  conditions  given  in  Tables  3.6-3.10  apply. 


f  <  1 

f  =  1 

f  >  1 

1. 

(200000000110) 

3 . 6a 

3 . 6a 

3.6b(f<l+/2) 

2. 

(020000000110) 

3 . 7a (f >0 . 5) 

3.7b 

3.7c 

3. 

(100001000110) 

3 . 8a 

3 . 8a,b 

3 . 8b  (f <2 ) 

4. 

(010001000110) 

3.9a  (f >0 . 5) 

3 . 9a , b 

3.9b 

5. 

(110000000110) 

3.10a  (f=0. 5) 

— 

3.10b  (f-l+/2) 
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Table  3.6a  0  £  f  ^  1  -2f(f  +1)  (200000000110) 

Adjacent  Vertex  Instability  Condition 


1.  (200000000011) 

f  x>2z 

2.  (100010000110) 

f  (f+i)  jA>(2f+D  jl 

3.  (100001000110) 

<f+i)jA>j2 

4.  (100000100110) 

2fV2j3 

5.  (100000010110) 

fjA>2j4 

6.  (100000001011) 

fjAX>V 

7.  (020000000110) 

2f (f+1)  <jA/jB)2  >(4f2  +  f=l) jA/jB+f (2f-l) 

8.  (002000000110) 

2f (f+1) (jA/jc)2>(4f2+3f+l) (jA/jc)2+f (2f+l) 

9.  (000200000110) 

f(f+l)(jA/jD-)2>jA/jD-+l 

10.  (010000001110) 

2(f+l)j2>(l-f)jBj5 

11.  (001000001110) 

2(f+l)j2>(l+2f)  jcj5 

12.  (000100001110) 

(f+1)^^D^5 

13.  (000010001110) 

'2(f+1^A^^5 

14.  (000001001110) 

2(f+1)j2>j2j5 

15.  (000000101110) 

2(f+1)j2>j3j5 

16.  (000000011110) 

(£+1)j!>2j4j5 

17.  (110000000101) 

f  (f+D  jAy  >2(i-f)  jBz 

18.  (101000000101) 

f (f+1) jAy >2(2f+l) jcz 

19.  (100100000101) 

f  (f+1)  jAy  >4  jD~z 

20.  (100010000101) 

f  (f+D  jAy  >23^ 

21.  (100001000101) 

f  (f+D  jAy  >2j2z 

22.  (100000100101) 

f  (f+D  jAy  >4  j  3z 

23.  (100000010101) 

f  (f+D  jAy  >8j4z 

24.  (010000001101) 

f (f+1) jAy  >(l-f ) jBj5z 

' 


. 


. 


Table  3.6a  (continued) 


25.  (001000001101) 

26.  (000100001101) 

27.  (000010001101) 

28.  (000001001101) 

29.  (000000101101) 

30.  (000000011101) 


f (f+l)j^y>(l+2f) jcj 
f(f+l)j2>2jD-j5z 
f  (f +  1 )  j  ?  >  j  ,  j  c-Z 

rl  X  J 

f (f+l)j^y>j2j5z 
f (f+1) j^y> j  3 j  5z 
f (f+1) j^y>4j4j5z 


« 
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Table  3.6b  f  >  1 

2  (f  -2f-l )  (200000000110) 

Adjacent  Vertex 

Instability  Condition 

1. 

(200000000101) 

| f 2-2f-l | y>2 (f-1) z 

2. 

(200000000011) 

| f 2-2f-l | x>4  f z 

3. 

(020000000110) 

2 | f 2-2f-l | ( jA/jB) 2>2f (f-3) (jA/jB)+f 

4 

(002000000110) 

|f2-2f-l| (jA/jc)2>2f (l+3f) (jA/jc)+f (l+2f) 

5. 

(000200000110) 

|f2-2f-l| (jA/jD-)>2(2f-l) (jA/jD-)+l 

6. 

(100010000110) 

|f2-2f-l| jA>(l+2f)j1 

7. 

(100001000110) 

|f2-2f-l| jA> ( 2— f ) j  2 

8. 

(100000100110) 

|f2-2f-l| jft>fj3 

9. 

(100000010110) 

|f2-2f-l| jA>4£j4 

10. 

(100000001110) 

2 | f 2-2f-l | jA>f (f-1) j5 

Table  3.7a  0.5  <  f 

<  1  f  (l-2f )  (020000000110) 

Adjacent  Vertex 

Instability  Condition 

1. 

(020000000011) 

( 2f-l ) x>  2z 

2. 

(020000000101) 

f (2f-l) y>2  (1-f) z 

3. 

(200000000110) 

f  (2f-l)  (jB/jA)2>(4f2+f-l)  (jB/jA)+2f (f+1) 

4. 

(002000000110) 

( 2  f — 1 )  ( j  B/ j  c ) 2  >  2 (1+f )  (jB/jc)  +  (l+2f) 

5. 

(000200000110) 

f(2f-l)(jB/jD-)2>jB/jD-+2 

6. 

(010010000110) 

f (2f-l) jfi>(2f+l) jx 

7. 

(010001000110) 

fVj2 

8. 

(010000100110) 

(2f-l) j  B>  j  3 

9. 

(010000010110) 

<2f"l) j  B>  4  j  4 

(2f-l) jB> (1-f ) j5 


10.  (010000001110) 


« 

\ 
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Table  3.7b  f  =  1  f(l-2f)  (020000000110) 


1. 

Adjacent  Vertex 

(020000000011) 

Instability  Condition 

(2f-l) x>2z 

2. 

(010010000110) 

f  (2f-l) jB>(2f+l) j1 

3. 

(010001000110) 

f  (2f-l) jB>(2f-l) j2 

4. 

(010000100110) 

(2f-D  jB>  j3 

5. 

(010000010110) 

(2f“l) j  B>  4  j  4 

6 . 

(010000001011) 

(2f-l) jgx>2j5z 

7. 

(200000000110) 

f  (2f-l)  (jB/jA)2>(4f2  +  f-l)  (jB/jA)+2f  (f+1) 

8. 

(002000000110) 

(2f-l)  (jB/jc)2>2(l+f)  (jB/jc)  +  (l+2f) 

9. 

(000200000110) 

f(2f-l)(jB/jD-)2>jB/jD-+2 

10. 

(001000001110) 

(2f-X)  j2>(l+2f)  jcj5 

11. 

(000100001110) 

(2f-l)j2>2jD-j5 

12. 

(000010001110) 

(2f-i)j|>j1js 

13. 

(000001001110) 

(2f-l) j|>j2j5 

14. 

(000000101110) 

(2f-l)j2>j3j5 

15. 

(000000011110) 

(2f-l)j2>4j4j5 

16. 

(011000000101) 

f  (2f-l) jBy>2 (f+2) jQz 

17. 

(010100000101) 

f  (2f-l) jBy>4jD-z 

18. 

(010010000101) 

f  (2f-l) jBy>2j1z 

19. 

(010001000101) 

f  (2f-l) jBy>2j2z 

20. 

(010000100101) 

f  (2f-l) jBy>2(3-f) jBz 

21. 

(010000010101) 

f  (2f-l) jBy>8j4z 

22. 

(001000001101) 

f (2f-l) j2y>2(l+2f) jcj5z 

23. 

(000100001101) 

f (2f-l)j2y>4jD-j5z 

24. 

(000010001101) 

f  (2f-l) j2y>2j  j  z 

. 


- 


Table  3 . 7b  (continued) 


25.  (000001001101) 

26.  (000000101101) 

27.  (000000011101) 


f  (2f-l) jgy>2j2j5z 
f  (2f-l) jgy>2j3j5z 
f  (2f-l) jgY>8j4j5z 


Table  3.7c  f  >  1 

Adjacent  Vertex 
1.  (020000000011) 

2.  (020000000101) 

3.  (200000000110) 

4.  (002000000110) 

5.  (000200000110) 

6.  (010010000110) 

7.  (010001000110) 

8.  (010000100110) 

9.  (010000010110) 

10.  (010000001110) 


-f  (020000000110) 

Instability  Condition 
x>  2z 

fy>2 (f-1) z 

f (jB/jA>  2>2f (f-3) jB/jA+2 (l+2f-f 2) 
<jB/jc)2>4f  (jB/jc)+(l+2f) 
f(jB/jD-)2>(2f-l)jB/jd-+2 
f jB> (4f-l) jx 

fVj2 
'  V33 

jB>(f-1>j5 
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Table  3.8a  0  v<  f  s<  1  -2f  (100001000110) 


Adjacent  Vertex 

Instability  Condition 

1. 

(100001000011) 

f  x>  2z 

2. 

(100001000101) 

fy>  2z 

3. 

(010001000110) 

2fja>±<l-2f)jBC£®J 

4. 

(001001000110) 

2fjA>(l+4f) jc 

5. 

(000101000110) 

fW 

6. 

(000011000110) 

fjA>2jl 

7. 

(200000000110) 

j2>(f+1)jA 

8. 

(000001100110) 

2fV*3 

9. 

(000001010110) 

fV2^4 

10. 

(000001001110) 

2^A^5 

Table  3.8b  1  ^  f  ^  2 

2  (f-2 )  (100001000110) 

Adjacent  Vertex 

Instability  Condition 

1. 

(100001000011) 

(2-f ) x>2z 

2. 

(100001000101) 

(2-f ) y>2z 

3. 

(010001000110) 

2<2-f)jA>jB 

4. 

(001001000110) 

2(2-f)jA>(l+4f)jc 

5. 

(000101000110) 

(2‘f>W 

6. 

(000011000110) 

(2-f) jA>2jx 

7. 

(200000000110) 

(2-f) j2> (l+2f-f2) jA 

8. 

(000001100110) 

2'2-f>jA>j3 

9. 

(000001010110) 

<2_f)jA>2j4 

10. 

(000001001110) 

2 (2-f) j A>  f  j 5 

' 
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Table  3.9a  0.5  <  f  ^  1  (l-2f)  (010001000110) 


1. 

Adjacent  Vertex 

(010001000011) 

Instability  Condition 

( 2f-l ) x>  2z 

2. 

(010001000101) 

( 2f-l ) y>  2z 

3. 

(100001000110) 

<2f-l)jB>2fjA 

4. 

(001001000110) 

(2f-l) jB>(l+4f) jc 

5. 

(000101000110) 

(2f-1)jB>2jD- 

6. 

(000011000110) 

(2f-l ) jB>  4 j 1 

7. 

(020000000110) 

j2>fjB 

8. 

(000001100110) 

(2f-l) jB> j3 

9. 

(000001010110) 

(2f-l)jB>4j4 

10. 

(000001001110) 

<2 f  1  >  j B>f  j  5 

Table  3.9b  f  >,  1 

-1  (010001000110) 

1. 

Adjacent  Vertex 

(010001000011) 

Instability  Condition 

x>2z 

2. 

(010001000101) 

y>2z 

3. 

(100001000110) 

jB>2(2-f)jA 

4. 

(001001000110) 

jB>(l+4f)jc 

5. 

(000101000110) 

jB>2V 

6. 

(000011000110) 

V4jl 

7. 

(020000000110) 

h>fH 

8. 

(000001100110) 

jB>j3 

9. 

(000001010110) 

2B>4^4 

10.  (000001001110) 


- 
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Table  3.10  a  f  =  0.5  Vertex  (110000000110)  ( l-f-4 f 2 ) 


Adjacent  Vertex 

Instability  Condition 

1. 

(200000000110) 

|  l-f-4f 2  |  jB>2f  (f+DjA 

2. 

(020000000101) 

| l-f-4f 2 | jAy>2 (1-f ) jB 

3. 

(020000000011) 

| l-f-4f2 I jAx>2f jBz 

4. 

(100001000110) 

| l-f-4f 2  I jB>2f j2 

5. 

(011000000110) 

| l-f-4f 2 | jA>2f (f+1) jc 

6 . 

(010100000110) 

U-f-«2|jA>jD- 

7. 

(010010000110) 

|l-f-4f2| jA>jx 

8. 

(010000100110) 

[l-f-4f2 | jA>f j3 

9. 

(010000010110) 

| l-f-4f 2 1 jA>4f j4 

10. 

(010000001110) 

| l-f-4f 2 | jA>f (l-f)j5 

Table  3.10b  f  =  1  +  /2~ 

Vertex  (110000000110)  2f 

Adjacent  Vertex 

Instability  Condition 

1. 

- (020000000110) 

2<f-3>jA>jB 

2. 

(200000000101) 

f  (3-f ) jBy>2 (f-1) jAz 

3. 

(200000000011) 

(3-f ) jBx>4jAz 

4. 

(101000000110) 

(3-f) jfi>(3f+l) jc 

5. 

(100100000110) 

f  (3-f ) jB> (f-1) jD- 

6. 

(100010000110) 

f  (3-f) jB> (l+2f) j1 

7. 

+ (100001000110) 

f  (3-f) jB> (f-2) j2 

8. 

(100000100110) 

(3-f)  j  b>:^  3 

9. 

(100000010110) 

(3-f) DB>4j4 

10. 

(100000001110) 

2(3-f) jB> (f-1) j5 

128. 


Secondly,  a  parallel  exists  between  each  of  the 
mixed  destabilizing  terms  and  the  corresponding  pure 
destabilizing  term  from  which  it  was  derived.  In  the 
region  f  >  1,  for  example,  if  one  replaces  one  of  the 
j^'s  in  an  adjacent  positive  vertex  of  the  pure  negative 
vertex  by  j2,  then  an  adjacent  positive  vertex  of  the 
mixed  negative  vertex  is  obtained.  The  same  situation 
holds  for  the  pure  B  destabilizing  term  and  the  mixed 
destabilizing  term  derived  from  it  when  f  >  1.  In  the 
region  f  >  1,  the  instability  conditions  for  each  pure 
destabilizing  term  are  identical,  in  all  cases  but  one, 
to  those  of  the  corresponding  mixed  destabilizing  terms, 
except  that  the  coefficient  \c+/C~\  is  modified.  The 
exception  occurs  in  the  instability  condition  which 
involves  the  currents  of  the  mixed  destabilizing  terms, 
(i.e.,  jA  and  j2  or  j  and  j2).  ln  this  case,  the 
inequality  involving  those  currents  in  the  pure  destabil¬ 
izing  terms  is  reversed  in  the  corresponding  mixed 

destabilizing  term  inequality.  For  example,  the  NEEC  B 
.  2 

pure  term  DBxy  has  the  instability  condition  fj  >  j2  in 
the  region  f  >  1,  and  the  corresponding  mixed  term  j^j^xy 
has  the  instability  condition  j0  >  fj  .  This  instability 

Z  13 

condition  implies  that  the  cones  of  negative  a 2  due  to 
each  of  these  destabilizing  terms  have  a  common  hyperplane 
and  that  if  fjg  >  j2  then  the  cone  due  to  the  pure 
destabilizing  term  applies.  Since  both  the  adjacent 
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vertices  are  negative  for  this  hyperplane,  it  does  not 
form  part  of  the  boundary  of  stability.  Furthermore,  any 
hyperplane  corresponding  to  the  remaining  instability 

.  .  o 

conditions  for  the  pure  term  j_,xy  lies  parallel  to  the 

a 

corresponding  hyperplane  in  the  mixed  term  j^joXy,  and 

J3  Z 

the  most  destabilizing  hyperplane  is  determined  by  the 

coefficient  |C  /C+ | .  This  coefficient  is  always  smaller 

for  the  pure  destabilizing  term  instability  conditions 

and  so  the  hyperplanes  for  this  term  destabilize  a 

greater  volume  of  parameter  space,  since  they  are  less 

2 

restrictive.  Similarly  the  pure  destabilizing  term  j  xy 

has  instability  conditions  which  destabilize  a  greater 

volume  of  parameter  space  than  the  mixed  term  j ^ j  2 XY - 

In  conclusion,  in  the  region  f  >  1  where  II  is  maximal, 

a2 

the  instability  conditions  due  to  a  pure  destabilizing 
term  in  determines  a  cone  of  negative  which  contains 
the  corresponding  cone  due  to  the  mixed  destabilizing 
term.  It  follows  that  when  the  current  in  the  EEC  2  is 
large  ( j 2  >>  jA  or  jfi)  then  the  cone  of  negative  a ^ 
defined  by  the  corresponding  pure  destabilizing  term  is 
replaced  by  the  cone  defined  by  the  mixed  destabilizing 
term,  and  so  the  unstable  region  in  parameter  space 
shrinks . 

In  the  region  f  <  1,  the  situation  is  the  same  for 

the  pure  and  mixed  destabilizing  terms  involving  the 

2 

NEEC  B.  However,  the  pure  destabilizing  term  j  xy  lies 
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in  a  truncated  region  of  II  when  f  <  1.  The  truncation 

a2 

of  II  due  to  the  critical  1-cycle  at  x  eliminates  the 
2 

positive  vertices  (200000000101)  and  (100000001110)  which 

are  adjacent  to  the  vertex  corresponding  to  the  pure  A 

destabilizing  term  for  f  >  1.  As  discussed  before,  this 

occurs  because  neither  the  NEEC  A  nor  the  EEC  5  have  an 

x-stabilizer ,  and  no  2-cycle  can  be  formed  between  the 

2 

parameters  of  the  terms  j^xz  and  3Aj5xY*  The  net  result 
is  that  the  instability  conditions  2j^  >  j,_  and  fy  >  2z 
present  in  Table  3.8a  for  the  mixed  destabilization  term 
jA32xY  cannot  be  formed  in  the  pure  A  destabilizing  term. 
It  follows  that  many  new  instability  conditions  appear 
for  the  pure  destabilizing  term  in  Table  3.6a.  However, 
if  the  missing  inequalities  given  above  are  multiplied 
by  the  remaining  instability  conditions  of  Table  3.8a, 
then  the  Table  3.6a  is  obtained,  apart  from  the  correct 
coefficients.  The  coefficients  \C+/C~\  of  the  instability 
conditions  due  to  the  pure  destabilizing  term  are  smaller 
than  the  coefficients  of  the  corresponding  mixed  destabil¬ 
izing  term,  and  hence  destabilize  a  greater  portion  of 
parameter  space.  It  follows  once  again  that  the  cone  of 
negative  due  to  the  mixed  destabilizing  term  is 
contained  in  the  cone  due  to  the  pure  destabilizing  term. 
This  result  holds  for  f  <  1 ,  in  the  case  of  the  NEEC  A. 
When  f  =  1  the  NEEC  B  also  has  a  critical  1-cycle  of  x, 
although  since  the  NEEC  A  equals  the  NEEC  B  at  f  =  1, 
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this  reduces  to  the  case  described  above  for  the  pure  A 
destabilizing  term  in  the  region  f  <  1. 

In  summary,  the  cones  of  negative  a2  due  to  the  mixed 
destabilizing  terms  are  contained  in  the  cones  defined 
by  the  corresponding  pure  destabilizing  terms.  When  the 
current  in  the  EEC  2  is  large,  the  cones  due  to  the  mixed 
destabilizing  terms  apply,  and  stabilize  the  unstable 
region  in  a2  as  defined  by  the  pure  destabilizing  terms. 

The  cone  of  negative  can  be  given  as  a  first  approxima¬ 
tion  by  the  union  of  the  cones  defined  by  the  pure 
destabilizing  terms.  In  this  case  the  instability 
condition  defined  between  the  currents  of  the  mixed 
destabilizing  term  is  dropped  from  the  complete  set  given 
for  the  corresponding  pure  destabilizing  term.  This 
approximation  would  overestimate  the  unstable  region  in 
a2  when  the  current  in  the  EEC  2  is  large,  but  can  be  used 
to  give  an  upper  bound  on  the  unstable  volume  of  parameter 
space  resulting  from  the  a2  instability. 


* 
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3.3.2.  The  Cone  of  Negative  a ^ 

When  f  ^  1 ,  all  of  the  terms  in  are  positive.  In 
the  region  f  >  1,  two  mixed  destabilizing  terms  occur  as 
follows : 

D  jgj3xyz  -2f ( f-1 ) 

2 )  jBj2j3xyz  -2  (f-1) 

The  diagrams  which  sum  to  form  these  terms  are  given  in 
figure  3.10,  along  with  the  coefficients  in  each  region 
of  f,  as  found  from  the  interaction  diagrams  in  figure 
3.8.  These  terms  have  the  positive  coefficients  6f  (1-f ) 
and  6  (i-f )  in  the  region  f  <  1,  for  the  first  and  second 
terms  above,  respectively. 

The  change  in  sign  of  these  terms  varies  in  a  regular 
way  with  the  parameter  f.  The  stabilizers  are  dominated 
by  the  inhibitory  positive  feedback  (PFI)  2-cycle  between 
x  and  y,  and  the  inhibitory  positive  feedback  (PFI) 

3-cycle  when  f  >  1.  At  f  =  1,  these  diagrams  have  coef¬ 
ficients  which  cancel  each  other.  When  f  <  1,  the  stabil¬ 
izers  are  greater  than  the  destabilizing  2  and  3-cycles. 

If  the  PFI  3-cycle  was  absent  these  terms  would 
be  positive  for  all  f  showing  that  it  is  the  addition  of 
this  3-cycle  to  the  PFI  2-cycle  which  causes  the  instab¬ 
ility  in  a^.  The  NEEC  A  also  forms  these  cycles  with 
the  EEC's  2  and  3,  but  has  a  positive  coefficient  for  all 
f.  Only  the  NEEC  B  in  the  region  f  >  1  forms  a  PFI 


Vertex  1  (020000100111) 

•  •  •  • 

a„  —  +  - 

J  •  • 

f  >  1  2f  (f-1)  +  2f  -  2f 

f  <  1  2f(l-f)  +  2f  -  2f 

Vertex  2  (010001100111) 

a3=  ••+••+•• 

0  0  0 

f  >  1  2f  +  2  (f-1)  +  2  - 

f  <  1  2f  +  2  (1-f )  +  2  - 


(j^xyz) 


-  2f 2  +  2f  =  -2f (f-1) 

-  2f 2  +  2f  -  6f  (1-f) 


^B^2^3Xyz) 


2f  -  2 f  ~  2  f  +  2  =  -2  (f-1) 

« 

2f  -  2f  -  2f  +  2  =  6  (1-f) 


Figure  3.10 

The  diagrams  contributing  to  the  II  negative  vertices 

a3 

in  the  region  f  >  1  for  the  reversible  model.  The 


coefficients  of  these  terms  are  positive  when  f  <  1. 


134  . 


2-cycle  and  PFI  3-cycle  which  are  large  enough  to  dominate 
the  positive  diagrams.  From  figure  3.10  it  can  be  seen 
that  these  destabilizing  cycles  increase  a  power  of  f 
faster  than  the  positive  diagrams  only  when  f  >  1.  It 
is  interesting  to  note  that  extreme  current  B  is  the 
only  NEEC  with  the  reverse  of  reaction  (4)  (except  B  , 

C  ,  D  which  were  eliminated  from  the  analysis).  This 
reaction  in  the  forward  direction  increases  the  value 
of  the  x-stabilizer  with  second  order  kinetics. 

One  of  the  most  significant  points  about  the 

destabilizing  terms  is  that  the  presence  of  the  EEC  3  is 

necessary  before  the  PFI  3-cycle  can  be  formed  with  the 

NEEC  B.  In  particular,  the  corresponding  pure  NEEC  term 
•3  .  .  .  . 

j^xyz  is  positive.  This  is  quite  different  from  the  a ^ 
case  discussed  previously,  in  which  thd  mixed  destabiliz¬ 
ing  terms  did  not  occur  without  a  corresponding  pure 
destablilizing  term. 

The  instability  conditions  which  define  the  cone 
of  negative  for  f  >  1  are  given  in  Table  3.11.  Several 
of  the  instability  conditions  in  Table  3.11  occur  in  the 
(*2  instability  conditions  for  the  destabilizing  terms 
involving  the  NEEC  B,  only  the  coefficients  are  modified. 
The  remaining  instability  conditions  are  new.  Therefore 
the  cone  of  negative  partly  overlaps  the  cone  of  nega¬ 
tive  o^.  It  follows  that  both  the  unstable  regions 
defined  by  and  are  necessary  to  describe  the 
boundary  of  the  unstable  steady  states  in  the  reversible 


model . 
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Table  3.11  The  instability  conditions  for  a 


3* 


(f  >  1) 

-2f  (f-1 )  (020000100111) 

Adjacent  Vertex 

Instability  Condition 

1. 

(200000100111) 

f (f-l) j^>2 (f 2-l) j2+4f (f+1) jAjB 

2. 

(100000101111) 

f(f-l)j2>(f2— l)jAj5 

3. 

(030000000111) 

VjB 

4. 

(020000001111) 

^3^5 

5. 

(010010100111) 

f (f-l) jB>(l+5f) jx 

6. 

(010001100111) 

7. 

(010000110111) 

(f-l)jB>4j4 

8. 

(002000100111) 

(f-l) (jB/jc) 2>4 (1+f ) jB/jc+3 (1+f ) 

9. 

(001000101111) 

(f-l) j|>3 (l+£) jcj5 

10. 

(000200100111) 

f (f-l)  (jB/jD-)2>f (jB/jD-)+2 

11. 

(000100101111) 

f  (f-l)  jB>  (1  +  f )  jD~j5 

12. 

(000010101111) 

f (f“l) jg>2 (1+f) j  X j  5 

13. 

(000000111111) 

(f'1)jB>434j5 

(f  >  1) 

-2  (f-l)  (010001100111) 

1. 

(200001000111) 

3b33>8^ 

2. 

(200000100111) 

(f-DjBj2>2(f2-l)j2 

3. 

(100001001111) 

jBV4jAj5 

4. 

(100000101111) 

<f-l) jBj2> (f2-l) jAj5 

5. 

(020001000111) 

6. 

(020000100111) 

)2>f)B 

7. 

(010001001111) 

j3>j5 

8. 

(001001100111) 

(f-l) jB>3 (1+f) jc 

9. 

(001000101111) 

(f-l) jBj2>3f (1+f) jcj5 

< 
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10. 

(000101100111) 

(f-l)jB>2jD- 

11. 

(000100101111) 

V2>(1+f)  vj 

12. 

(000011100111) 

(f-1) jB>4j1 

13. 

(000010101111) 

(f-1) jBj2>2(l+f) jxj 

14. 

(000001110111) 

15. 

(000000111111) 

<f-l) jBj2>4f j4 j5 
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3.3.3.  The  Cone  of  Negative 

The  destabilizing  terms  in  a ^  generate  many  negative 

vertices  in  II.  through  the  product  a,a».  However,  as 

discussed  in  section  3.2.2,  it  is  the  NMO  negative 

vertices  which  have  a  negative  feedback  (NF)  3-cycle 

passing  through  a  critical  1-cycle  that  may  extend  the 

unstable  region  in  the  network  by  an  appreciable  amount. 

In  the  region  f  >  1  there  are  eight  negative  vertices 

in  II.  ,  and  none  of  these  are  NMO  vertices.  The  exponent 
2 

polytope  nA  has  maximum  configuration  in  this  region  of 

f,  except  at  f  =  2  and  1  +  /2.  At  f  =  1  +  /2  the 

32  32  22  22 

destabilizing  terms  jAx  y,  jAxy  ,  jAj2x  Y  and  ^2xy 

.  2 

vanish  in  because  the  coefficient  of  the  a 2  term  jAxY 

is  zero.  The  polytope  II ^  is  truncated  and  two  new  nega- 

2  .  2  2 
tive  vertices  appear,  corresponding  to  the  terms  3A3BX  Y 

2  2 

and  j  j 3xy  ,  which  are  on  an  edge  of  the  maximum  polytope. 

These  terms  have  MOD'S  and  do  not  destabilize  the  network 

beyond  the  cone  of  negative  •  At  f  =  2,  no  new  negative 

vertices  appear,  in  parallel  with  a  results.  The 

negative  vertices  in  IT.  are  listed  in  figure  3.11  for 

2 

f  >  1. 

In  the  region  f  <  1,  there  are  seventeen  negative 

vertices  in  II.  .  This  large  increase  in  the  number  of 

2 

negative  vertices  occurs  because  the  critical  1-cycle  in 

the  NEEC  A  causes  many  of  the  MO  vertices  in  II.  to 

2 

vanish.  At  the  special  case  of  f  =  0.5,  the  coefficient 
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2 

of  the  a0  term  j  xy  vanishes  and  the  polytope  II.  is 
ZB  A2 

further  truncated  because  the  MO  vertices  corresponding 

32  32  22  22 

to  the  terms  jgx  y,  jfixy  ,  jgj2x  y  and  jgj2xy  vanish. 

No  new  destabilization  occurs  in  the  network,  however, 

•  2  2 

because  the  new  II.  vertices  j,j„x  y  and  j^j^xy  which 

A9  Aj B  2  JAJB  2 

^  3  3 

are  negative  lie  on  the  j  -j  edge  in  n.  which  lies 

A  B  A2 

inside  the  region  of  instability  determined  by  a  2 ■  The 
other  new  vertices  which  appear  are  all  positive  as  the 
a.  term  jDx  still  exists  (i.e.,  the  x-stabilizer  does 

_L  ID 

3 

not  vanish) .  The  NMO  term  j_.xyz  is  positive  and  corresponds 

ID 

to  a  midpoint  of  an  edge  of  II.  at  f  =  0 . 5 ,  as  expected. 

z 

When  f  =  1,  the  critical  1-cycle  in  the  NEEC  B 

vanishes  and  the  polytope  II.  is  truncated  further  still. 

z 

For  this  value  of  f,  24  negative  vertices  exist.  However, 
this  is  artificial  since  current  A  equals  current  B  and 

i 

the  polytope  II.  is  symmetric  about  the  hyperplane 

a2 

•  3  3 

through  the  3A“DB  edge. 

The  appearance  of  new  negative  vertices  in  II.  with 

2 

f  is  shown  in  figure  3.11,  although  the  special  cases 
f  =  0.5,  2  and  1  +  /2  are  omitted.  Only  when  the  poly¬ 
tope  is  truncated  from  its  maximum  configuration  at  f  >  1 
due  to  a  critical  1-cycle,  do  NMO  negative  vertices  appear. 
Although  the  polytope  is  also  truncated  for  f  =  0.5, 

2  and  1  +  /2  when  the  coefficient  of  an  a ^  destabilizing 
term  vanishes,  no  additional  NMO  negative  vertices  are 

formed.  The  truncation  of  II.  that  occurs  for  these 

A2 


. 


. 


.  > 

• 
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values  of  f  is  not  destabilizing  since  the  stabilizers 
present  before  the  truncation  all  exist  and  contribute 
to  the  diagrams  which  define  the  positive  adjacent 
vertices  to  the  negative  vertices.  The  new  negative 
vertices  that  did  appear  were  midpoints  on  a  destabiliz¬ 
ing  edge  before  truncation  occurred.  Hence  they  cannot 
define  hyperplanes  which  lie  outside  the  unstable  region 
of  the  maximal  polytope. 

The  negative  vertices  which  occur  in  II.  are  listed 

z 

in  Table  3.12.  The  coefficients  in  this  table  were 
calculated  by  summing  all  the  possible  diagrams  which 
can  be  formed  in  the  product  The  possible 

diagrams  contributing  to  this  product  are  given  in  the 
interaction  diagrams  in  figure  3.8.  Only  the  product 

is  needed  to  calculate  the  MO  vertices,  while  the 
NMO  vertices  have  diagrams  that  must  be  subtracted 
from  the  diagonal  product  a^c^-  0nly  when  the  negative 
feedback  3-cycle  in  survived  total  cancellation  did 
negative  NMO  vertices  result. 

Two  negative  NMO  vertices  occur  in  II.  throughout 

2 

the  region  f  <  1,  and  they  correspond  to  the  following 
destabilizing  terms: 

j^xyz  -8f  (1+f )  (Vertex  3,  table  3.12) 

^A-^5X^Z  -4f  (1+f )  (Vertex  23,  table  3.12) 

Both  of  these  terms  are  due  to  a  negative  feedback  (NF) 


c 


i 


*  i  -  ■ - -  M 


\ 
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Table  3.12  Negative  Vertices 


term 

1.  (300000000210) 

2.  (300000000120) 

3.  (300000000111) 

4.  (210000000210) 

5.  (030000000210) 

6.  (030000000120) 

7.  (030000000111) 

8.  (201000000210) 

9.  (200100000210) 

10.  (200010000210) 

11.  (200001000210) 

12.  (021000000210) 

13.  (020100000210) 

14.  (020010000210) 

15.  (020001000210) 

16.  (100002000210) 

17.  (100002000120) 

18.  (010002000210) 

19.  (010002000120) 

20.  (200000100210) 

21.  (200000010210) 

22.  (200000001120) 

23.  (200000001111) 

24.  (020000100210) 

25.  (020000010210) 

26.  (020000001120) 

27.  (020000001111) 


coefficient  f  <  1 
0 

-2  (1+f ) 2f 
-8f  (1+f) 

-2f  (1+f)  (1-f) 
f  (2f2-3f+l) 
f2  (l-2f ) 

-2f  ( 3 f — 2 ) 

-2f  (2f+l )  (f+1) 
-4f  (1+f ) 

-2f  (f+1) 

-2f  (f+1) 

-3f  (2f 2-l ) 

-4  f 2+f +1 

-4f 2+2f +1 

l-2f 

-  2f 

-2f 

(l-2f) 

( 1— 2f ) 

-2f  (f+1) 

-8f  (f+1) 

-2f 2  (f+1) 

-4f  (f+1) 
f  (2-3f ) 

4f  (2-3f ) 
f2 ( 2  —  3  f ) 

-2f  ( 5f-4 ) 


f  >  1 

2f 3-6f 2+2f +2 
4f  (f 2-2f-l ) 

16f  (f-2 ) 

2  (2f2-7f2  +  4f+l) 
f (1-f) 

-f2 

-2f  (2-f ) 
9f3-9f2-9f-3 
2  (4f2-7f-l) 

2  (3f2-3f-2) 

2 (2f2-5f+l) 
f (4f2-6f-l) 

2f 2-5f +1 
4f 2-2f +1 
l-2f 
2  (f-2) 

-4 

-1 

-1 

5f 2-8f-l 
8  (2f2-3f-l) 

2f  (2f2-3f-l) 

8f  (2f-3 ) 
f (f-2) 

4f  (f-2) 
f 2  (f-2 ) 

-2 f  ( 4  —  3 f ) 


■ 

* 

< 


7 

- 

-*>  <  ,  :  A  )  '  • 


Table  3.12  (continued)  Sign  in  the  Region  of  Dominance. 


f =0 

0<f <0 . 5 

f =0 . 5  0 . 5<f <1 

f=l 

2>f  >1 

f  >2 

2.414 

1. 

0 

0 

0  0 

0 

— 

- 

+ 

2. 

0 

- 

—  — 

— 

— 

- 

+ 

3. 

0 

- 

—  — 

— 

4. 

0 

— 

_  — 

0 

5. 

0 

+ 

_(D 

0 

— 

— 

— 

6. 

0 

+ 

_  (1  2 ) 

— 

— 

— 

— 

7. 

0 

— 

8. 

0 

- 

—  — 

— 

9. 

0 

— 

—  - 

— 

10. 

0 

— 

—  — 

- 

11. 

0 

— 

—  — 

- 

12. 

0 

— 

13. 

+1 

— 

14. 

+1 

— 

15. 

+1 

— 

16. 

0 

- 

_  - 

— 

— 

+ 

+ 

17. 

0 

— 

—  — 

— 

- 

— 

- 

18. 

+1 

0 

— 

— 

— 

— 

- 

19. 

+1 

0 

— 

— 

— 

— 

— 

20. 

0 

— 

- 

— 

21. 

0 

— 

—  — 

- 

22. 

+4 

- 

-  - 

— 

23. 

0 

— 

—  — 

— 

24. 

+1 

— 

25. 

0 

— 

26. 

0 

— 

27. 

0 

— 

(1)  At  f  =  0.5  vertex  5  is  replaced  by  (120000000210) 

(2)  At  f  =  0.5  vertex  6  is  replaced  by  (120000000120) 
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3-cycle  which  passes  through  a  critical  1-cycle.  The 
diagrammatic  simplification  for  these  terms  is  shown  in 
figure  3.12.  The  critical  1-cycle  at  x  in  the  NEEC  A 
(f  <  1)  prevents  the  cancellation  of  the  NF  3-cycle, 
which  originates  in  a^,  by  the  product  a^t^.  In  t*ie  mixed 
current  vertex,  the  NEEC  A  (f  <  1)  overlaps  with  the 
EEC  5  which  is  the  only  other  extreme  current  without  an 
x-stabilizer .  Again,  this  prevents  the  formation  of  a 
stabilizer  diagram  in  the  product  . 

The  instability  conditions  which  define  the  cone  of 

negative  A2  for  the  NMO  negative  vertices  above  are  given 

in  Table  3.13.  The  hyperplanes  defined  by  the  adjacent 

positive  vertices  are  identical  to  hyperplanes  found  in 

2 

the  instability  conditions  for  the  a ^  term  j  xy  when 
f  <  1  (Table  3.6a).  The  remaining  instability  conditions 
in  Table  3.13  are  due  to  adjacent  negative  vertices  with 
MOD's  and  do  not  form  part  of  the  boundary  of  stability. 

As  in  the  case,  the  instability  conditions  for  the 
A2  NMO  terms  can  be  multiplied  to  give  the  additional 
inequalities  which  appear  in  Table  3.6a.  The  cones  of 
negative  defined  by  the  NMO  vertices  3  and  23  do  not 
extend  the  cone  of  negative  defined  by  the  destabiliz- 
ing  term  j  xy.  When  the  current  m  the  EEC  2  is  large 
and  the  boundary  of  stability  in  a  2  is  described  by  the 
mixed  destabilizing  term  j^j2xy,  the  network  cannot  be 
further  destabilized  by  either,  since  the  cones  of 


* 
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Vertex  3  (300000000111)  j  3xyz 

A 2  =  a!a2  -  a3 


Vertex  23  (200000001111)  j  ^j_xyz 

A  J 


_ ^ 


2  =  + 


-4f  (1+f ) ,  f  1 


Vertex  7  (030000000111)  j_3xyz 

13 


•  • 


_  ^ 


A2  - 


■2f(3f-2),  f=l 


Vertex  27  (020000001111) 


jB2j5xyz 


-2f(5f-4),  f=l 


Figure  3.12 

The  diagrams  contributing  to  the  n.  NMO  negative  vertices 

A2 

in  the  region  f  <  1  for  the  reversible  model. 


. 

-  . 


Table  3.13 


The  instability  conditions  for  the  NMO 
vertices  in  the  reversible  model. 


A2 

NMO  Vertex  3 

Adjacent  Vertex 

(300000000111) 

Coefficient 

-  8f(f  +  1)  f  <  1 

Instability  condition 

1. 

(300000000012) 

16  (1+f ) 

f  z>  2z 

2. 

(300000000120) 

-2  (l+f)2f 

4z> (1+f ) y 

3. 

(210000000210) 

-2f  (1+f)  (1-f ) 

4zjA> (1-f )xjB 

4. 

(201000000210) 

-2f  (2f+l )  (1+f) 

4z j A> (2f+l)xjc 

5. 

(200100000210) 

-4f  (1+f ) 

2zjA>XjD" 

6. 

(200010000210) 

-2f  (1+f) 

4zVxji 

7. 

(200001000210) 

-2f  (1+f) 

4zjA>xj2 

8. 

(200000100210) 

-2f  (1+f) 

4zVX^3 

9. 

(200000010210) 

-8f  (1+f) 

ZVX*4  ' 

10. 

(210000000120) 

16  (1-f) 

f (f+l)yjA>2 (1-f ) zjB 

11. 

(201000000102) 

16 (2f+l) 

f (f+l)yjA>2(2f+l)zj 

12. 

(200100000102) 

32 

f (f+l)yjA>4zjD- 

13. 

(200010000102) 

16 

f (f+l)yjA>2zj1 

14. 

(200001000102) 

16 

f (f+l)yjA>2zj2 

15. 

(200000100102) 

32 

f (f+l)yjA>4zj3 

16. 

(200000010102) 

64 

f (f+l)yjA>8zj4 

17. 

(200000001111) 

-4f  (1+f) 

2jA>j5 

*' 
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Table  3.13  (continued) 


&2  NMO  Vertex  23 

(200000001111) 

-  4f  (f+1) 

1.  (300000000111) 

-8f  (1+f) 

j5>2jA 

2.  (210000000210) 

-2f  (1+f)  (1-f) 

2z j  5>  (1-f )xjB 

3.  (201000000210) 

-2f  (2f +1 )  (1+f) 

2z j  5> (2f+l)xjc 

4.  (200100000210) 

-4f  (1+f) 

ZVXV 

5.  (200010000210) 

-2f  (1+f) 

2z j  5>x3 i 

6.  (200001000210) 

-2f  (1+f) 

2z j  5>xj  2 

7.  (200000100210) 

-2f  (1+f) 

2z j  5>xj  3 

8.  (200000010210) 

-8f  (1+f) 

z^5>xj4 

9.  (200000001120) 

-2f 2  (1+f) 

2z>fy 

10.  (100000002012) 

4 (1+f) 

fxVzj5 

11.  (010000002102) 

4 (1-f) 

f (1+f ) y j A> (1-f) zjBl5 

12.  (001000002102) 

4 (l+2f ) 

f  (1+f  )yj  +  (l+2f )  zjcj 

13.  (000100002102) 

8 

f (l+f)yj2>2zjD-j5 

14.  (000010002102) 

4 

f(i+f)yjl>zj1j5 

15.  (000001002102) 

4 

f (l+f)yj^>zj2j5 

16.  (000000102102) 

4 

f (l+f) y j^>z j  3 j  5 

17.  (000000012102) 

16 

f  d+f)yj^>4zj4j5 

' 

. 


. 


■ 

negative  defined  in  Table  3.13  lie  inside  the  cone 

2 

due  to  the  pure  destabilizing  term  j  xy. 

When  f  =  1  a  critical  1-cycle  in  the  NEEC  B  appears, 
allowing  the  NF  3-cycle  present  in  this  extreme  current 
to  dominate  in  Two  new  negative  NMO  vertices  appear 

corresponding  to  the  following  destabilizing  terms: 


^BXyZ 

-2 

(Vertex 

7,  table  3.12) 

3Bj5XyZ 

-2 

(Vertex 

27,  table  3.12) 

The  diagrams  corresponding  to  these  NMO  terms  are  also 

.  3 

given  in  figure  3.12.  Note  that  the  NMO  term  jgXyz  is 

the  same  NMO  destabilizing  term  which  occurred  at  f  =  1 

in  the  irreversible  model.  For  f  <  1  the  stabilizer 

diagrams  from  ouo^  a re  only  partially  cancelled,  and  the 

terms  become  positive  for  f  <  2/3  and  f  <  4/5.  However, 

these  terms  are  no  longer  vertices  of  II.  when  f  <  1. 

2 

Since  the  extreme  currents  A  and  B  are  equal  at  f  =  1. 
these  NMO  vertices  need  not  be  considered  since  they  are 
equivalent  to  the  negative  NMO  vertices  involving  current 
A  previously  described. 


■ 

/ 

.  v> 
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3.3.4.  Summary 

The  region  of  instability  in  the  reversible  model  is 
a  union  of  the  cones  of  negative  a2  and  for  f  >  1, 
and  is  given  by  a  union  of  the  cones  of  negative  a 2  and 
f°r  f  <  1.  The  dominant  cycles  of  the  destabilizing 
terms  are  summarized  in  Table  3.14,  along  with  the  range 
of  f  in  which  these  cycles  destabilize  the  network  and 
the  extreme  currents  involved.  The  unstable  range  of 
f  is  f  >  0  in  the  reversible  model. 

In  the  reversible  model,  there  are  two  pure  destabil¬ 
izing  terms  in  a2  which  arise  from  the  NEEC ' s  A  and  B. 

In  addition  there  are  two  mixed  destabilizing  terms 
resulting  from  the  overlap  of  current  A  or  B  with  the 
EEC  2.  The  cone  of  negative  due  to  either  of  the  pure 
destabilizing  terms  contains  the  cone  due  to  the  correspond¬ 
ing  mixed  destabilizing  term.  It  follows  that  when  the 
current  in  the  EEC  2  is  larger  than  that  in  the  NEEC ' s  A 
and  B  by  a  certain  amount,  the  unstable  region  in  a2  is 
determined  by  the  mixed  destabilizing  terms.  This  reduces 
the  unstable  volume  of  parameter  space  that  would  be 
found  if  the  pure  destabilizing  terms  were  used  as  an 
approximation  to  the  unstable  region  of  a 2  alone.  The 
two  NMO  destabilizing  terms  in  A2  present  at  f  <  1  do 
not  extend  the  boundary  of  the  unstable  steady  states  as 
defined  by  ct2,  and  so  are  not  necessary  to  describe  the 
unstable  region  in  the  network. 
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Table  3.14 

The  destabilizing  cycles  found  in  the  reversible  model 
analysis . 


Destabilizing  cycles  Extreme  currents 

involved 


Range  of  f 
(unstable) 


_  ^ 


/ 


A2  (NMO) 


A 


-  -  ^  * 


pure  A 

f 

< 

1  +  /2 

pure  B 

f 

> 

0.5 

mixed 

A, 

2 

f 

< 

2 

mixed 

B, 

2 

f 

> 

in 

• 

o 

mixed 

B  / 

3 

f 

> 

1 

mixed 

B , 

3,  2 

f 

> 

1 

pure  A 

f 

1 

pure  B 

- 

f 

= 

1 

mixed 

A, 

5 

f 

1 

mixed 

B , 

5 

f 

— 

1 

* 


A  critical  1-cycle  at  x  was  necessary  in  addition. 


. 


. 


In  the  region  f  >  1,  two  mixed  destabilizing  terms 


appear  in  a^.  The  inhibitory  positive  feedback  3-cycle, 

which  is  necessary  for  the  destabilization,  requires 

the  presence  of  the  EEC's  to  be  formed.  It  follows  that 

the  destabilizing  terms  cannot  be  formed  without  the 

reverse  reactions  and  that  this  destabilization  will  be 

absent  in  the  irreversible  model,  as  the  results  of 

section  3.2.3  confirm.  C  extends  the  boundary  defined 

a3 

by  C  only  for  certain  hyperplanes  and  values  of  f  (>1) . 
a2 

This  extension  depends  on  the  value  of  the  coefficients 
of  the  destabilizing  terms  in  and  a^.  For  example, 
take  the  instability  condition  on  the  currents  and 
j^.  The  inequalities  are: 


2 *  ^ B  >  4j4 

(table 

\ 

3.7c) 

3 :  >  434 

(table 

3.11) 

In  the  region  1  v<  f  <  2,  C  determines  the  boundary; 

a2 

for  f  >  2  C  determines  the  boundary.  At  f  =  2,  C  =  C 
a3  a2 

on  this  hyperplane,  i.e.,  in  summary 


1  <  f  <  2 
f  =  2 
f  >  2 


C 

C 

c 


a 

a 

a 


2 

2 

2 


The  a2  instability  condition  is  less  restrictive  close  to 
f  =  1  and  therefore  destabilizes  a  greater  portion  of 
parameter  space  in  this  region  of  f. 


. 


' 


\ 
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The  destabilizing  cycles  found  in  this  analysis 
support  previous  findings,  are: 

o^:  inhibitory  positive  feedback  2-cycle 

2  8 

(competition  interaction  ) 

:  inhibitory  positive  feedback  2-cycle 

2  8 

inhibitory  positive  feedback  3-cycle 

2  8 

A2  (NMO) :  negative  feedback  3-cycle 

Field  and  Noyes^  found  the  same  unstable  range  of  f  in 

the  irreversible  model:  0.5  <  f  <  1  +  /2 .  The  extension 

of  unstable  f  in  the  reversible  model  was  later  calculated 
2  8 

by  Field  to  be  0  f  ^  1.51  for  the  rate  constants  of 

step  (5)  in  the  model  (1.6)  given  by  the  arbitrarily 

-1  -5  -1 

chosen  values  k^  =  1  sec  and  k_^  =  10  sec 
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3.4.  Re-evaluation  of  the  omitted  extreme  currents 

The  extreme  currents  A  ,  B  ,  C  and  D  were  eliminated 
from  the  preceeding  analysis  of  the  reversible  model,  to 
reduce  the  number  of  parameters  and  because  the  unstable 
region  was  thought  to  lie  close  to  the  AB  edge  of  the 
current  polytope  II  .  However,  on  the  basis  of  the  destab¬ 
ilization  found  in  a^,  it  is  possible  that  mixed  destabil¬ 
izing  terms  involving  the  EEC's  could  be  formed  with  the 
omitted  NEEC ' s  above.  In  ,  the  mixed  destabilizing 
terms  cannot  form  unless  a  corresponding  pure  destabiliz¬ 
ing  term  is  present.  Only  in  is  a  mixed  destabilizing 
term  without  a  corresponding  pure  destabilizing  term 
possible . 

The  interaction  diagrams  for  the  extreme  currents 
A  ,  B  ,  C  and  D  are  given  in  figure  3.13  for  each  region 
of  f.  The  only  destabilizing  cycle  present  is  the  PFI 
2-cycle  between  x  and  z.  As  expected,  the  pure  terms 
are  positive  in  the  region  f  >  1  and  positive  or  zero  in 
the  region  f  <  1.  Similarly  the  pure  terms  are 
positive  or  zero.  The  absence  of  a  critical  1-cycle 
prevents  the  destabilization  of  by  an  NMO  terms  involv¬ 
ing  the  negative  feedback  3-cycle  found  in  the  NEEC ' s  C 
and  A  (f  >  1) .  In  conclusion  no  pure  destabilizing 
terms  can  be  formed  in  the  NEEC ' s  A  ,  B  ,  C  and  D. 

The  results  of  the  reversible  model  analysis  suggest 
that  the  only  possible  source  of  destabilization  involving 
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Figure  3.13 

The  interaction  diagrams  of  the  extreme  currents  A  ,  B  ,  C 
and  D  in  the  reversible  model. 
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NEEC's  A-,  b“  and  C~  is  the  formation  of  mixed  inhibitory 

positive  feedback  (PFI)  3-cycles  in  a^.  The  EEC  2  has  a 

inhibitory  positive  feedback  (PFI)  2-cycle  which  can 

contribute  an  inhibitor  (dashed  arrow)  to  form  a  PFI 

3-cycle  with  any  of  these  NEEC's.  The  EEC  3  is  necessary 

in  addition  for  another  mixed  PFI  3-cycle  in  a^.  When 

.  2  . 

the  diagrams  contributing  to  each  of  the  terms  j^j2xYz 
and  j.j9j._xyz,  where  i  =  A  ,  B  or  C  ,  were  summed,  only 
the  NEEC's  B~  and  C-  in  the  region  f  <  1,  were  found  to 
result  in  a  negative  coefficient.  The  diagrams  correspond¬ 
ing  to  the  destabilizing  terras 


^B-j2XyZ 

2  (f-1) 

f 

< 

1 

V^3XyZ 

2  (f-1 ) 

f 

< 

1 

.2  . 

3c-32xyz 

2  (3f-l) 

f 

< 

1/3 

Vj2j3XyZ 

2  (3f-l ) 

f 

< 

1/3 

found  in  the  region  f  <  1  are  given  in  figure  3.14.  It 
is  interesting  to  note  that  both  NEEC's  B  and  C  have 

the  reverse  of  reaction  (4) .  The  NEEC  B  (f  >  1)  which 

* 

was  responsible  for  the  destabilizing  terms  in  section 
3.3.2  also  has  the  reverse  of  reaction  (4).  It  can  be 
seen  that  reaction  (4)  increases  the  value  of  the  x- 
stabilizer ,  and  so  has  a  stabilizing  effect,  when  present 
in  an  NEEC.  Both  NEEC's  B~  and  C~  in  the  region  f  <  1 
form  destabilizing  mixed  terms  in  which  positive  feedback 
(PFI)  2  and  3-cycles  are  independent  of  f,  while  the 
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2 

Destabilizing  term  j  -j_xyz 

B  Z. 


f<l=4+2f-4-2=  2  (f-1 ) 
f>l=4f+2f-4-2=  6  (f-1) 

Destabilizing  term  j j  2  j  3XY z 


f<  l=2+2f+4-4-2-2=  2  (f-1) 
f>l=2+2f+4f  -4-2-2=  6  (f-1) 

.  .  .  .  2  . 

Destabilizing  term  3 ^_j 2xyz 


=  2  (f +2 )  +2f-4  +  2f-2  =  2  (3f-l)  -  ,f  <  1/3 

Destabilizing  term  jc“j2j3xyz 


=  2  +  2  (f+2 )  +2f+2f  -2-4-2=  2(3f-l),  -,f<l/3 

Figure  3.14 

The  diagrams  contributing  to  the  mixed  destabilizing  terms 
in  for  f  <  1 ,  in  the  reversible  model,  corresponding  to 
the  omitted  NEEC ' s ,  B  ,  C  . 
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positive  diagrams  have  f  dependence.  This  allows  the 
PFI  2  and  3-cycles  to  dominate  when  f  <  1.  Note  also  the 
direction  of  the  PF  3-cycles  in  the  terms  involving  B 
and  C  are  reversed  in  direction  with  respect  to  the 
3-cycle  in  the  B  (f  >  1)  term. 

The  method  of  destabilization  found  in  in  the 
region  f  <  1  is  the  same  as  the  previously  found  case  in 
the  region  f  >  1.  There  is  one  difference,  however,  and 
that  is  the  lack  of  a2  destabilization  by  the  NEEC ' s  C  , 
B-.  In  fact,  the  NEEC's  C~  and  B~  do  not  form  any  pure 
destabilizing  terms  in  cl^,  a ^  or  This  means  that 

in  the  current  polytope  II for  the  reversible  model,  an 
unstable  region  appears  on  the  edge  joining  the  NEEC  B 
(f  <  1)  and  the  EEC  2  or  the  NEEC  C_  (f  <  1)  and  the 
EEC  2,  although  all  the  extreme  points  of  II  correspond¬ 
ing  to  these  extreme  currents  are  stable  steady  states. 
When  the  EEC  3  is  included,  an  isolated  volume  of  unstable 
steady  states  occurs  in  nc,  which  does  not  extend  to  any 
of  the  extreme  points. 

It  is  significant  that  the  EEC's  are  necessary  for 
this  instability,  and  that  there  is  no  relationship 
between  the  a3  destabilization  in  the  region  f  <  1  and 
the  irreversible  model.  In  fact,  if  all  the  reactions  in 
the  irreversible  model  are  simply  reversed,  the  Oregonator 
is  stable  for  all  f  values. 
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The  cones  of  negative  due  to  the  destabilizing 

terms  in  the  region  f  <  1  can  be  estimated  by  inspection 

of  the  instability  conditions  found  previously  for  the 

instability  in  the  region  f  >  1  (Table  3.11).  The 

resulting  sets  of  instability  conditions  are  shown  in 

Tables  3.15  a  to  d.  The  actual  cones  of  negative 

will  always  lie  inside  the  cones  defined  by  the  instability 

conditions  of  Tables  3.15  a  to  d,  because  if  an  adjacent 

positive  vertex  of  IT  has  been  omitted  the  actual  cone 

a3 

of  negative  is  smaller. 
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Table  3.15a  (020001000111)  2(f-l) 


0  f  <  1 


Adjacent  Vertex  Instability  Condition 


1. 

(002001000111) 

(1-f) j|> (l-3f ) j2 

2. 

(000201000111) 

3. 

(021000000111) 

(1-f) j2>f (f+3) jc 

4. 

(020100000111) 

<l-f)j2>4fjD 

5. 

(020010000111) 

(1-f) j2>(l+f) jx 

6. 

(020000010111) 

(1-f) j2>4fj4 

7. 

(010011000111) 

(1-f) jB>4j1 

8. 

(010001100111) 

V^3 

9. 

(010001010111) 

(l-f>V4j4 

10. 

(010001001111) 

V^s 

11. 

(110001000111) 

VjA 

Table  3.15b  (010001100111) 

Adjacent  Vertex 

2  (f-1)  0  4  £  <  1 

Instability  Condition 

1. 

(020001000111) 

'  j3>jB 

2. 

(020000100111) 

V 

3. 

(001001100111) 

(1-f) jB> (3f-l) jc 

4  o 

(000101100111) 

(l-f)jB>6jD 

5. 

(000011100111) 

(1-f) jB>4j1 

6. 

(000001110111) 

(1-f ) jB>4 d4 

7. 

(100000101111) 

2jBj2>(f+2)jAj5 

8. 

(100001001111) 

iB33>2jAj5 

9. 

(010001001111) 

j3>j5 

10. 

(001000101111) 

d-f )  jBj2>f  (1+f )  ^C^5 

11. 

(000100101111) 

(i-f)  jBj2>(l+5f)  jDj5 

12. 

(000010101111) 

(l-f)jBj2>4(l+f)j1j5 

13. 

(000000111111) 

d-f)  jBj2>4fj4j5 

■  ‘ 
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Table  3.15c  (002001000111)  2  (3f-l)  0  N<  f  <  1/3 


1. 

Adjacent  Vertex 

(020001000111) 

Instability  Condition 

<l-3f) j£>(l-f) j2 

2. 

(101001000111) 

VjA 

3. 

(000201000111) 

(l-3f) j|>3j^ 

4. 

(003000000111) 

(l-3f) j2>f (1+f) jc 

5. 

(001011000111) 

(l-3f) jc>4j1 

6. 

(001001100111) 

jC>j3 

7. 

(001001010111) 

(l“3f) jc>4 j4 

8. 

(001001001111) 

(l-3f) jc>(l+f) j5 

9. 

(201000000111) 

(l-3f ) jcj 2>f jA 

Table  3.15d  (001001100111) 

Adjacent  Vertex 

1.  (002001000111) 

2  (3f-l)  0  <:  f  <  1/3 

Instability  Condition 

j3>jC 

2. 

(002000100111) 

(l-3f) j2>f (1+f) jc 

3. 

(010001100111) 

(l-3f) jc>  (1-f ) jB 

4. 

(000101100111) 

(l-3f ) jc>6jD 

5. 

(000011100111) 

(l-3f) jc>4j1 

6. 

(000001110111) 

d-3f )  jc>4j4 

7. 

(100000101111) 

2(l-3f) jcj  2> (1-f)  (f+2) jAj 

8. 

(100001001111) 

(l-3f) jcj3>2 (X-f ) jAj5 

9. 

(001001001111) 

(l~3f ) j  3> (1+f) j  5 

10. 

(001000101111) 

(l-3f) j2>f (1+f) j5 
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3.5.  Comparison  of  the  reversible  and  irreversible  model 
results . 

The  effect  of  reversibility  in  a  network  is  demonstrated 
by  the  change  in  the  complete  set  of  extreme  currents  when 
the  reverse  reactions  are  added  to  the  irreversible  net¬ 
work.  All  the  extreme  currents  present  in  an  irreversible 
network  are  of  the  NEEC  type  and  appear  in  the  correspond¬ 
ing  reversible  network.  As  a  result,  an  instability  if 
present  in  an  irreversible  network  always  appears  in  the 
reversible  case  as  well.  If  new  instability  is  to  occur 
in  the  reversible  network,  it  must  involve  the  extreme 
currents  of  the  reversible  network  which  cannot  be 
formed  in  the  irreversible  case.  Any  such  extreme  current 
will  contain  at  least  one  reversed  reaction.  The  new 
extreme  currents  in  the  reversible  Oregonator  model 
consist  of  the  EEC's  and  the  new  NEEC's  which  were  formed 
by: 

i)  Reversing  all  the  reactions  in  one  of  the  NEEC's 
of  the  irreversible  model,  for  example:  A  (f  >  1) , 

B_  (f  <  1)  and  C  . 

ii)  Mixing  different  forward  and  reverse  reactions  to 
form  new  NEEC's,  for  example:  D,  D  ,  A  (f  <  1) 
and  B  (f  >  1)  . 

Note  that  a  reaction  and  its  reverse  can  only  appear 
together  in  an  EEC,  and  every  reaction  step  in  a  network 
forms  an  EEC  of  the  network. 


'< 
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The  results  of  the  analysis  performed  on  the 
Oregonator  model  are  most  easily  compared  at  f  =  1,  since 
in  this  case  the  extreme  currents  A  and  B  are  equal. 

The  instability  in  the  irreversible  model  at  f  =  1  is 
due  to  one  pure  destabilizing  term  in  a2  which  is  derived 
from  the  NEEC  B  (=A) .  In  the  reversible  model  the  same 
pure  destabilizing  term  occurs  in  at  f  =1,  while  the 
new  NEEC's  B~ ,  C~ ,  D,  D~  and  the  EEC's  1  to  5  cause  no 
additional  pure  destabilizing  terms.  As  discussed 
previously  the  pure  NMO  destabilizing  term  in  A2  does 
not  extend  the  unstable  region  in  either  model  beyond  the 
boundary  defined  by  a2.  The  new  mixed  destabilizing 
terms  in  a2  and  A2  of  the  reversible  model  also  define 
cones  of  instability  which  do  not  extend  the  unstable 
region  defined  by  the  pure  destabilizing  term  in  a2 . 

At  f  =  1  the  reversible  model  instability  conditions 
for  the  pure  and  mixed  destabilizing  terms  in  a2  contain 
the  irreversible  model  instability  conditions.  Therefore, 
the  a2  destabilizing  terms  in  the  reversible  model  at 
f  =  1  cannot  extend  the  unstable  region  defined  by  the 
irreversible  model  hyperplanes.  The  instability  condi¬ 
tions  for  the  reversible  model  will  define  hyperplanes 
between  the  new  current  parameters  in  addition  to  the 
hyperplanes  found  in  the  irreversible  model  instability 
conditions.  The  new  hyperplanes  of  the  reversible  model 
cause  further  bounds  on  the  unstable  region  of  parameter 


, 
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space.  For  example,  in  the  reduced  set  of  extreme 
currents  used  in  the  reversible  model  analysis,  the 
rate  of  the  reverse  reactions  (3)  and  (5)  are  given  by 
j  3  and  respectively.  Hence  the  new  instability 
conditions  involving  these  currents  such  as 
(Table  3.7b,  f  =  1)  are  equivalent  to  conditions  on  the 
rates  of  the  reverse  reactions  necessary  for  instability. 

When  f  differs  from  one  the  situation  changes, 
although  the  destabilizing  terms  in  both  cases  (reversible 
vs.  irreversible)  still  contain  the  NEEC ' s  A  or  B.  The 
destabilizing  extreme  currents  in  the  irreversible  model 
are  the  NEEC's  A  (f  >  1)  and  B  (f  <  1) .  If  the  reactions 

are  allowed  to  reverse,  then  NEEC  B  can  be  extended  into 

the  region  f  >  1,  and  NEEC  A  can  be  extended  into  the 

region  f  <  1,  by  reversing  reactions,  (4)  and  (1) 
respectively.  This  extension  of  the  NEEC's  A  and  B  in 
the  reversible  model  is  diagrammed  in  figure  3.15.  These 
new  NEEC's  cannot  be  combined  with  the  corresponding 
irreversible  model  NEEC's  since  a  reaction  when  reversed, 
alters  the  interaction  diagram.  Compare,  for  example, 
the  NEEC's  A  (f  <  1)  and  A  (f  >  1)  in  figure  3.15  with 
their  interaction  diagrams  in  figure  3.8.  In  general, 
when  an  arrow  has  a  coefficient  which  is  a  function  of 
a  varying  parameter,  such  as  f,  so  that  the  arrow  changes 
sign  at  some  value,  the  extreme  current  has  to  be  consid¬ 
ered  separately  in  each  region. 
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irreversible  model 


reversible  model 


Figure  3.15 

The  extension  of  the  NEEC's  A  and  B  in  the  reversible  model. 
The  heavy  line  indicates  the  reaction  which  is  reversed  in 
the  irreversible  model  extreme  current,  when  the  new  destab¬ 
ilizing  extreme  current  is  formed  in  the  reversible  model. 


The  new  NEEC's  A  (f  <  1)  and  B  (f  >  1)  each  define 
a  pure  destabilizing  term  in  cl^  of  the  reversible  model. 

As  a  result  of  these  terms,  the  unstable  range  of  f  in 
the  reversible  model  is  extended  below  0.5  and  above 
1  +  /2,  the  limits  for  unstable  behavior  in  the  irrevers¬ 
ible  model.  The  question  now  arises  as  to  whether  these 
terms  can  also  extend  the  bounding  hyperplanes  of  the 
unstable  region  for  f  within  the  unstable  range  of  the 
irreversible  model.  To  answer  this  question,  consider 

the  three  irreversible  model  instability  conditions  in 

.2 

a.2  corresponding  to  the  pure  destabilizing  terms  j^xy(f  >  1 
and  j^xy(f  <  1).  These  instability  conditions  are 
present  in  the  reversible  model  also,  however  the  coef¬ 
ficients  for  the  same  hyperplanes  are  modified  for  the 
destabilizing  terms  j^xy(f  <  1),  j^xy(f  >  1),  j^2xY  an<^ 
j_,j0xy.  Table  3.16  lists  the  three  irreversible  model 
instability  conditions  and  the  coefficient  as  a  function 
of  f,  as  they  appear  in  the  four  destabilizing  terms  in 

(*2  of  the  reversible  model. 

When  a  plot  of  \C+/C~\  versus  f  is  made  as  in  figure 
3.16,  the  term  in  a2  which  determines  the  most  destabiliz¬ 
ing  hyperplane  can  be  found.  The  coefficient  of  each 
destabilizing  term  corresponding  to  one  of  the  hyperplanes 
labelled  (1),  (2)  and  (3)  in  Table  3.16  is  plotted  in 

figure  3.16a,  b  and  c  respectively.  The  label  R  means 
that  the  curve  corresponds  to  a  coefficient  which  occurs 


■ 
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Table  3.16 

Comparison  of  the  irreversible  model  hyperplanes  found  in  the 
four  a 2  destabilizing  terms  of  the  reversible  model. 

f  <  1  |c+/c~| 


(1)  j2xy  (Table  3.6a) 

1.  l/2f (1+f ) (jA/jc)2>l/2(l+3f+4£2) jA/jc+f (l+2f) 

2. 

3.  fx>2z  2/f 


(2)  jAj2xy  (Table  3.8a) 

1.  fjA>(l+4f)jc 

2.  fy>2z 

3.  fx>2z 

(3)  j._.xy  (Table  3.7a)  (irreversible) 

1.  (2f-l ) (jB/jc)2>(l+2f)+2(l+f)jB/jc 

2.  f (2f-l)y>2  (l-f)z 

3.  (2f-l) x>2z 

(4)  jgj  2xy  (Table  3.9a) 

1.  (2f-l) jB> (l+4f) jc 

2.  (2f-l)y>2z 

3.  (2f-l ) x>2z 


l+4f/f 

2/f 

2/f 


2  ( 1-f ) /f ( 2f-l ) 
2/  (2f-l) 

l+4f/(2f-l) 

2/ (2f-l) 

2/  (2f-l) 


(continued  on  next  page) 
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Table  3.16  (continued) 


f  >  1  1 

C+/C" 1 

(1)  j2xy  (Table  3.6b)  (irreversible) 

1.  1/2 | f 2-2f-l | (jA/jc)2>f d+2f)+f (1+ef) jA/jc 

2.  |f2-2f-l|y>2(f-l)z 

2  (f-1 ) 

| f 2-2f-l | 

3.  | f2-2f-l|x>4fz 

4f 

1 f 2-2f-l 1 

(2)  2xy  (Table  3.8b) 


1.  (2-f) jA> (l+4f) jc 

(l+4f) 

(2-f) 

2.  (2-f )y>2z 

2/2-f 

3.  (2-f) x>2z 

2/2-f 

(3)  jBxy  (Table  3.7c) 

1.  (jB/jc)2>4f  (jB/jc)  +  (l+2f) 

2.  fy> 2  (f-1 ) z 

2 (f-l)/f 

3.  x>2z 

2 

(4)  jBj2xY  (Table  3.9b) 


i.  d+4f ) jc 

l+4f 

2.  y>2z 

2 

3.  x>2z 

2 

2 


. 

. 


167. 


Plot  of  the  coefficients  in  Table  3.16  versus  f ,  where  (a), 
(b) ,  (c)  correspond  to  the  first,  second  and  third  hyper¬ 

planes  in  Table  3.16,  respectively. 
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in  the  reversible  model  only,  while  an  IR  means  that  the 
coefficient  is  the  same  as  in  the  irreversible  model 
hyperplane .  If  the  coefficient  of  a  hyperplane  is  smaller 
than  the  coefficient  of  another  hyperplane  with  the  same 
parameters,  then  the  hyperplane  with  the  smaller  coef¬ 
ficient  is  less  restrictive  and  hence  destabilizes  a 
greater  portion  of  parameter  space.  The  parameters  x, 
y,  and  z  are  the  reciprocal  steady  state  concentrations 

X  ,  Y  ,  Z  and  so  if  x  >  y  then  Y  >  .  However,  the 

o  o  o  oo 

relative  relationship  between  parallel  hyperplanes  is 
the  same,  since  all  the  inequalities  are  inverted  in 
(x  ,  Yo,  Zq)  space.  Therefore,  it  is  sufficient  to 
determine  the  most  destabilizing  hyperplane  in  the  x,  y, 
z  parameters. 

The  first  thing  to  note  is  that  the  NEEC's  A  and  B  are 
equal  at  f  =  1,  and  so  no  extension  of  the  unstable  region 
defined  by  the  irreversible  model  can  occur  due  to  the 
pure  a2  destabilizing  terms  of  the  reversible  model. 
Secondly,  the  mixed  destabilizing  terms  in  the  reversible 
model  always  lie  either  on  or  inside  the  boundary  defined 
by  their  corresponding  pure  destabilizing  term.  The 
only  extension  of  the  unstable  region  defined  by  the 
irreversible  model  that  can  take  place  in  the  reversible 
model  due  to  a2,  is  related  to  the  pure  terms  involving 
the  new  NEEC's  A(f  <  1)  and  B(f  >  1).  It  follows  that 
the  mixed  destabilizing  terms  present  at  f  =  1  cannot 
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destabilize  the  reversible  model  and  so  no  extension 
occurs  at  f  =  1  for  any  hyperplane,  as  stated  in  the 
beginning  of  this  section. 

If  the  instability  conditions  in  Table  3.16  which 
are  independent  of  the  current  parameters  are  considered, 
several  comparisons  can  be  made.  In  the  region  f  <  1, 
the  unstable  range  of  f  is  extended  below  0.5  in  the 
reversible  model.  The  second  instability  conditions  for 
the  destabilizing  term  j^xy  (f  <  1)  in  Table  3.16  (f  <  1) 
is  missing,  but  other  inequalities  in  Table  3.6a  yield 
restrictions  on  the  ratio  of  y/z ,  which  will  either  be 
the  same  or  less  restrictive  than  the  corresponding  mixed 
destabilizing  term  instability  condition.  It  can  be 
seen  from  the  coefficient  of  the  hyperplane  (R2)  in 
figure  3.16b  that  this  hyperplane  is, more  destabilizing 
than  the  corresponding  irreversible  model  hyperplane  (R3) 
for  f  <  0.67.  In  other  words,  close  to  f  =  1  the 
instability  condition  found  in  the  irreversible  model 
determines  the  most  destabilizing  hyperplane.  The  third 
instability  condition  in  Table  3.16  for  f  <  1  due  to  the 
NEEC  A  (f  <  1)  defines  a  hyperplane  (Rl)  in  figure  3.16c, 
which  is  less  restrictive  than  the  corresponding  hyper¬ 
plane  (IR3 )  in  the  irreversible  model,  for  all  f  <  1. 
Therefore,  the  unstable  region  of  parameter  space  is 
extended  by  this  hyperplane  in  the  reversible  model,  within 
the  unstable  range  of  f  for  the  irreversible  model,  as 
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well  as  below  it.  However,  as  f  -*■  0,  the  stabilizing 
effect  of  the  positive  dominant  terms  which  occur  in  the 
remaining  reversible  instability  conditions  increases 
rapidly,  because  the  EEC's  define  terms  which  are 
independent  of  f,  while  the  NEEC  A  (f  <  1)  has  f^  dependenc 
(EEC  5  has  f  dependence  only.)  From  figure  3.16c  it  can 
be  seen  that  the  hyperplane  (Rl)  is  very  restrictive  for 
f  <  0.5,  since  the  coefficient  increases  rapidly  as  f  0. 
Hence,  the  unstable  volume  of  parameter  space  shrinks 
as  f  0. 

In  the  region  f  >  1,  the  new  NEEC  B(f  >  1)  causes 
the  unstable  range  of  f  present  in  the  irreversible 
model  to  be  extended  into  the  region  f  >  1  +  /2.  Two  of 
the  instability  conditions  in  Table  3.16  are  independent 
of  the  current  parameters.  The  instability  condition 
(3)  defines  a  hyperplane  (R3)  in  figure  3.16c,  involving 
the  NEEC  B ( f  >  1),  which  has  a  smaller  coefficient  and 
so  is  less  restrictive  than  the  corresponding  hyperplanes 
defined  by  any  of  the  other  destabilizing  terms.  In 
this  case  the  new  hyperplane  in  the  reversible  model 
destabilizes  a  greater  portion  of  parameter  space  for 
f  >  1,  and  the  bounding  hyperplane  of  the  unstable  region 
in  the  irreversible  model  is  extended.  The  second 
instability  condition  due  to  NEEC  B(f  >  1)  defines  the 
hyperplane  (R3)  in  figure  3.6b  which  extends  the  unstable 
region  beyond  the  corresponding  irreversible  model  hyper¬ 
plane  ( IRl )  only  for  f  >  1.6.  Therefore,  there  is  no 
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extension  of  the  unstable  region  defined  by  this  hyperplane 
in  the  irreversible  model  for  f  close  to  one. 

When  the  first  instability  condition  of  each  destabil¬ 
izing  term  in  Table  3.16  are  compared,  the  trend  found 
for  the  concentration  parameters  is  continued.  By 
referring  to  figure  3.16a,  it  can  be  seen  that  the  hyper¬ 
plane  defined  by  the  pure  destabilizing  terms  destabilize 
a  greater  portion  of  parameter  space  than  the  correspond¬ 
ing  mixed  destabilizing  terms  and  that  the  hyperplanes 
defined  by  the  terms  involving  the  new  NEEC  current 
parameters,  A(f  <  1)  and  B(f  >  1),  determine  the  most 
destabilizing  hyperplane  in  the  regions  f  <  1  and  f  >  1, 
respectively.  It  is  also  evident  that  the  coefficients 
of  the  most  destabilizing  terms  increase  rapidly  as 
f  ->  0  and  f  -*-  °°,  shrinking  the  unstable  volume  of  para¬ 
meter  space. 

In  summary,  the  new  NEEC ' s  A(f  <  1)  and  B(f  >  1) 
have  introduced  instability  into  not  present  in  the 

irreversible  model.  In  addition  to  extending  certain  of 
the  bounding  hyperplanes  defined  by  the  irreversible 
model  in  parameter  space,  the  range  of  unstable  f  itself 
was  increased  from  0.5  <  f  <  1  +  /2  to  f  >0.  However, 
the  region  of  instability  in  parameter  space  rapidly 
diminishes  as  f  drops  below  f  =  0.5  due  to  the  NEEC  A(f  <  1 
destabilizing  term  increasing  as  f  .  Also,  as  f  increases 
above  f  =  1  +  /2,  stabilization  of  the  destabilizing 
NEEC  B (f  >  1)  term  again  occurs,  since  the  adjacent 
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NEEC  A (f  >  1)  dominant  term  is  positive  for  f  >  1  +  /2 
and  increases  as  f  2 .  In  addition  the  NEEC  C  increases 
as  f2  so  that  the  destabilization  hyperplane  (R3)  in 
figure  3.16a  becomes  highly  restrictive,  and  only  for  a 
very  large  current  in  the  NEEC  B(f  >  1),  will  the  network 
be  unstable. 

As  discussed  in  section  3.3.4,  the  mixed 
destabilizing  terms  in  the  reversible  model  are  also 
necessary  for  describing  the  unstable  region  of  the 
network.  The  a3  instability  requires  the  presence  of 
EEC’s,  and  so  does  not  occur  in  the  irreversible  model 
for  any  value  of  f.  The  a3  instability  conditions  involve 
only  the  current  parameters  and  the  only  inequality 
that  can  be  compared  with  the  ones  in  Table  3.16  is 

(f-1) jB/jc)2>4 (1+f )  (jfi/jc)  +  3  (1+f )  (Table  3.11) 

which  involves  the  new  NEEC  B(f  >  1).  The  corresponding 
hyperplane  for  this  instability  condition  does  not 
destabilize  the  a2  hyperplane  (R3)  in  figure  3.16a 
until  f  >  1  +  /2.  For  example,  at  f  =  2.5,  j B/ j c  >  10.03 
while  the  coefficient  at  f  =  1.1  is  84.73.  The  large 
f  values  required  for  the  a3  instability  to  be  important 
follow  from  the  factor  (1-f)  contained  in  the  coefficient 
of  this  destabilizing  term.  Note  that  the  equivalent 
factor  (f-1)  appears  in  the  NEEC  B(f  <  1)  a3  destabilizing 
terms  discussed  in  section  3.4.  Since  the  experimentally 
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25 

determined  range  of  f  is  0  <  f  <  2,  the  instability 
present  in  the  region  f  >  1  is  uninteresting  experimentally. 
In  addition,  the  NEEC ' s  involved  in  the  a3  instability 

all  contain  the  reverse  of  reaction  (4),  which  has  a  very 

,,  „  ,  -10  -2  -1  28. 
small  rate  constant  (1.6  x  10  M  sec  .  ) 

So  far,  the  discussion  in  this  section  has  dealt 
with  the  shift  in  the  hyperplanes  of  the  irreversible 
model  when  the  reverse  reactions  are  included  in  the 
network.  However,  this  change  in  the  boundary  of  the 
unstable  steady  states  is  small  in  comparison  to  the 
possible  effect  of  new  hyperplanes  defined  by  the  remain¬ 
ing  reversible  model  instability  conditions,  when  f  is 
within  the  unstable  range  of  the  irreversible  model. 

These  new  hyperplanes  involve  the  new  extreme  currents 
which  have  a  stabilizing  effect  on  the  network  since 
they  place  additional  requirements  on  the  unstable  NEEC 
current  parameters  that  must  be  fulfilled  for  a  steady 
state  to  be  unstable.  These  requirements  are  missing 
in  the  irreversible  network,  and  so  the  validity  of  the 
irreversible  model  when  0.5  <  f  <  1  +  /2,  as  an  approx¬ 
imation  to  the  reversible  model  rests  on  the  restrictiveness 
of  the  new  hyperplanes.  There  are  new  hyperplanes 
involving  the  currents  of  each  of  the  EEC's  for  every 
destabilizing  term  in  the  reversible  model.  If  the  current 
in  either  of  the  EEC's  1,  3,  4  or  5  (f  >  1)  is  large 
(large  Keq.)  then  many  of  the  unstable  steady  states 
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accessible  in  the  irreversible  model  would  be  unattainable 
in  the  reversible  network,  since  the  current  in  the 
unstable  NEEC ' s  A  and  B  would  have  to  be  very  large  in 
order  to  destabilize  the  network.  This  condition  need 
not  be  satisfied  by  the  EEC  2.  As  discussed  in  section 
3.3.1,  if  the  current  in  the  EEC  2  is  large  then  the 
reversible  network  is  destabilized  by  the  cones  of  nega¬ 
tive  a2  due  to  the  mixed  destabilizing  terms  in  and 

although  these  cones  lie  inside  the  cones  defined  by  the 
corresponding  pure  destabilizing  a2  terms,  the  boundary 
of  the  unstable  region  is  not  severely  restricted. 
Similarly,  for  f  <  1  a  large  current  in  the  EEC  5  will 
still  allow  the  destabilization  of  the  network  by  the 
mixed  destabilizing  term  in  A2  involving  the  EEC  5.  In 
conclusion  the  irreversible  model  is  ,a  good  approximation 
of  the  reversible  model  when  the  currents  in  the  EEC's 
which  are  not  involved  in  any  mixed  destabilizing  terms 
in  the  network  are  small,  and  f  is  within  0.5<f<l+/?. 

A  summary  of  the  destabilizing  extreme  currents  in 
the  reversible  and  irreversible  models  is  given  in  Table 
3.17,  together  with  the  corresponding  range  of  f. 
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Table  3.17 

Summary  of  the  destabilizing  extreme  currents  in  the 
irreversible  and  reversible  models. 


Extreme  irreversible 
Current 

Reversible 

Destabilizing 

Cycle 

A  (f > 1 )  l£f<l+/2 

Uf<l+/2 

PFI  2-cycle 

B  (f <1 )  0.5<f£l 

0 . 5<f-$l 

“2 

A  (f  <  1) 

0<fSl 

B  ( f  >  1 ) 

f£l 

A  (f >1)  ,2 

l£f<2 

B  (f  >1)  /  2 

f>l 

no  extra 

A  (f <1)  ,  2 

0<f£l  1 

instability 

B  ( f  <  1 )  ,  2 

0. 5<f^l 

B  (f >1)  ,  3 

f>l 

PFI  2-cycle 

B(f>l) ,2,3 

f>l 

PFI  3-cycle 

a-3 

B~  (f <1 )  ,2 

0<f<l 

3 

B~  (f <1 )  ,2,3 

0<f<l 

C“(f<l)  ,2 

0-$f  <1/3 

C"(f<l) ,2,3 

0-^f  <1/3 

A=B  f=l 

f=l 

|  NF  3-cycle 

A2 

A  (f <1) 

0<f  ^1 

with  a 
critical 

A  (f <1)  , 5 

0<f^l 

1-cycle  no 
'  extra 

instability 
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3.6.  General  Conclusions 

In  any  chemical  system  there  are  only  so  many  ways 
of  taking  a  linearly  independent  sum  of  the  reactions 
with  positive  coefficients,  which  preserves  the  steady 
state  concentrations  of  the  internal  species.  Each  of 
these  ways  corresponds  to  an  extreme  current  of  the 
system,  and  any  general  steady  state  of  the  system  is  a 
positive  linear  combination  of  the  extreme  currents. 

Two  broad  classes  of  extreme  currents  can  be  identified 
in  a  reversible  network:  the  equilibrium  extreme 
currents  (EEC's)  and  the  non-equilibrium  extreme  currents 
(NEEC's).  In  the  preceding  analysis,  the  specific 
relationships  between  instability  and  the  extreme  currents 
were  discussed  for  the  reversible.  Oregonator  model.  The 
more  general  consequences  of  the  analysis  are  now 
summarized  by  considering  the  destabilizing  cycles  that 
were  necessary  for  a  cone  of  instability  in  the  network. 

The  pure  destabilizing  terms  which  occur  in  a 
reversible  network  are  derived  from  the  NEEC's  and  the 
destabilizing  cycles  are  the  same  as  those  given  in 
section  2.5  as  a  result  of  previous  work  on  irreversible 
networks.  The  fact  that  the  unstable  range  of  the 
stoichiometric  parameter  f  was  extended  in  the  reversible 
model  analysis  by  new  NEEC's,  points  out  that  an  irrevers¬ 
ible  network  can  be  completely  stable  and  yet  the  network 
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is  unstable  when  the  reactions  are  allowed  to  reverse. 

The  discussion  of  the  previous  section  pointed  out  that 
the  unstable  volume  of  parameter  space  is  small  in  such 
cases ,  and  that  the  current  in  the  responsible  unstable 
NEEC  must  be  large  for  an  unstable  steady  state  to  form. 

The  mixed  destabilizing  terms  which  appear  in  a 
reversible  network  all  have  a  current  component  from  a 
NEEC,  which  may  or  may  not  form  a  pure  destabilizing 
term  alone.  The  mixed  destabilizing  terms  in  a2  do  not 
occur  without  a  corresponding  pure  destabilizing  term 
derivable  from  the  same  NEEC.  The  rule  for  the  forma¬ 
tion  of  such  a  mixed  destabilizing  term  is  as  follows: 

A  network  which  is  destabilized  by  an  inhibitory 
positive  feedback  (PFI)  2-cycle  in  a2  will  form  a 
new  mixed  destabilizing  term  whose  dominant  cycle 
is  a  mixed  current  PFI  2-cycle,  when  an  EEC  with 
an  overlapping  PFI  2-cycle  is  present  in  the 
network . 

Such  a  mixed  destabilizing  term  stabilizes  the  cone  of 
instability  defined  by  the  corresponding  pure  destabiliz¬ 
ing  term  when  the  current  in  the  EEC  is  larger  than  the 
current  in  the  NEEC  by  a  value  determined  by  the  ratio 
of  the  coefficients  of  the  pure  and  mixed  destabilizing 
terms.  If  the  irreversible  model  of  a  network  has  an 
instability  in  a2  due  to  a  PFI  2-cycle,  then  the  cone  of 
instability  produced  by  the  same  PFI  2-cycle  in  the 
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reversible  network  will  be  comparable  as  long  as  the 
currents  in  the  EEC's  which  do  not  have  the  same  PFI 
2-cycle  are  small. 

The  situation  differs  in  k  >  2,  where  a  mixed 

destabilizing  term  may  exist  without  a  corresponding 
pure  destabilizing  term  being  possible.  The  character¬ 
istic  mixed  destabilizing  cycle  in  such  terms  is  a  mixed 
PFI  k-cycle,  k  >  2,  which  occurs  in  addition  to  a  PFI 
2-cycle.  This  combination  of  destabilizing  cycles  can 
be  dominant  in  a^,  k  >  2.  The  required  EEC's  for  a 
mixed  destabilizing  term  in  are  those  which  can 
contribute  an  inhibitor  to  the  PFI  k-cycle.  This  form 
of  destabilization  can  be  summarized  as  follows: 

Given  a  reversible  model  which  contains  a  NEEC 
with  a  PFI  2-cycle,  which  may  or  may  not  destabil¬ 
ize  the  network  in  pure  form,  then  the  formation 
of  mixed  destabilizing  terms  in  a^,  k  >  2,  is 
possible,  if  the  only  k-cycle  which  forms  is  a 
mixed  PFI  k-cycle  with  one  or  more  inhibitors 
contributed  by  the  EEC's. 

This  form  of  destabilization  occurs  only  in  reversible 
networks.  If  the  NEEC  involved  has  a  pure  stabilizing 
term  in  a,  ,  because  it  can  form  a  negative  feedback 

X. 

(NF)  k-cycle,  then  the  cone  of  instability  due  to  the 
mixed  current  PFI  k-cycle  will  be  important  only  if  the 
currents  in  the  EEC's  involved  in  the  PFI  k-cycle  are 
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large.  For  example,  in  the  reversible  Oregonator  model, 

.  2  • 

the  mixed  destabilizing  terms  for  f  >  1  are 

and  j  j  j  _xyz .  The  NEEC  B(f  >  1)  has  a  NF  3-cycle  which 

stabilizes  the  pure  B  term  in  for  f  >  1: 


xyz 


When  the  current  in  the  EEC  3  is  large  enough,  the  mixed 
PFI  3-cycle 


2 


/ 

B 


j^xyz 


or 


3Bj2j3xyz 


will  dominate  and  the  steady  state  becomes  unstable. 

The  inhibitor  in  the  EEC  3  which  causes  the  PFI  3-cycle 
above  is  a  result  of  the  reverse  of  reaction  (3) .  The 
corresponding  instability  conditions  found  in  are: 


VjB 


2 

for  the  destabilizing  term  jfij3xyz  and 


2>3 

3>:l 


in  the  case  of  j B j  2 □ 3XYZ • 

As  a  final  case,  the  formation  of  mixed  destabilizing 
terms  in  A2  will  be  discussed.  This  rule  is  actually  an 
extension  of  the  theorem  given  in  section  2.5,  which 
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gave  the  conditions  under  which  a  pure  destabilizing 
term  can  be  formed  in  A2°  The  formation  of  mixed  destab¬ 
ilizing  terms  in  A 2  does  not  occur  without  a  corresponding 
pure  destabilizing  term  as  in  the  a.^  case,  and  can  be 
summarized  as  follows: 

If  a  network  is  destabilized  in  A2  by  a  NF  3-cycle 

which  passes  through  a  critical  1-cycle,^'*'  then  an 

EEC,  which  has  a  vanishing  stabilizor  at  the  same 

parameter  as  the  critical  1-cycle,  produces  a  new 

mixed  destabilizing  term  which  has  a  mixed  NF 

3-cycle  which  differs  only  by  the  substitution  of 

one  arrow  from  the  original  pure  destabilizing  cycle. 

This  rule  will  also  hold  in  the  case  that  NF  k-cycles  in 

A,  n  destabilize  a  network.  The  cone  of  instability 
k-1 

determined  by  the  mixed  destabilizing  term  in  A2  will  be 
important  in  the  reversible  model  when  the  current  in 
the  EEC  involved  is  large. 
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